Visualisation des résultats : tableaux et graphiques
Laurie Auclair, Jocelyn Lauzon & Sabrina Brisson
10/09/2023
Importer les données
library(datasets)
data(iris)Exploration de données
Avant de débuter des analyses, il est important de se familisariser avec son jeu de données afin d’avoir une idée de sa structure. Cette étape permet d’identifier des motifs, des tendances ou des relations préalablement aux tests statistiques, puis de vérifier la qualité des données. Elle permet aussi d’évaluer si les données ont besoin d’être transformées avant de procéder aux analyses.
Quelques fonctions de base
Pour vérifier le nombre de lignes et de colonnes de votre dataset:
dim(iris)## [1] 150 5Pour voir les 6 premières lignes du jeu de données:
head(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosaPour connaître le nom des colonnes:
names(iris)## [1] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width" "Species"Pour vérifier la classe de chaque variable dans la table de données:
str(iris)## 'data.frame': 150 obs. of 5 variables:
## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...Pour obtenir certaines statistiques descriptives de base, telles que le minimum, le maximum puis la moyenne (variables continues), ainsi que le nombre d’observations (variables catégoriques):
summary(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100
## 1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300
## Median :5.800 Median :3.000 Median :4.350 Median :1.300
## Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199
## 3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800
## Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500
## Species
## setosa :50
## versicolor:50
## virginica :50
##
##
## Pour oconnaître les niveaux d’une variable catégorique:
levels(iris$Species)## [1] "setosa" "versicolor" "virginica"Visualisation d’un jeu de données
La fonction plot() est une fonction versatile qui permet de créer une grande variété de figures avec les données brutes. Dépendemment du type de donnée et des arguments fournis à la fonction, plot() peut générer plusieurs types de figures comme:
Une seule matrice de scatterplots figure peut être créée avec toutes les variables continues.
plot(iris)
Pour créer un scatterplot de la relation entre une variable spécifique et une autre, il suffit d’entrer celles-ci dans la fonction, séparées par le signe ~. Par exemple, si on veut visualiser la relation entre la longueur des sépales et celle des pétales dans le jeu de données iris, on peut l’écrire ainsi:
plot(iris$Sepal.Length~iris$Petal.Length)
À l’aide de la fonction boxplot(), on peut visualiser la dispersion d’une variable selon dles groupes d’une variable catégorique. Ce type de figure permet aussi d’identifier rapidement des valeurs aberrantes ou des anomalies.
boxplot(iris$Petal.Length~iris$Species)
La fonction hist() génère un histogramme de la variable qui y est précisée.
hist(iris$Sepal.Width)
Conseils pour des figures et tableaux réussis
Les figures sont des représentations visuelles des résultats. Elles rendent la lectures des résultats principaux plus facile et permettent de mettre en évidence des tendances ou motifs intéressants. On veut pouvoir la comprendre sans avoir à osciller entre la figure et le texte. Chaque figure devrait correspondre à un message général de votre étude. Chacune devrait répondre en totalité ou en partie à un objectif.
Voici quelques aspects à considérer pour présenter vos figures:
Un titre descriptif: doit contenir les variables mesurées, les unités de mesure, le nom commun et en latin du taxon (si applicable). Le titre doit fournir assez d’information pour qu’on comprenne le contexte de la figure sans devoir se référer au texte.
Titre des axes: doivent comprendre les variables et leur unité de mesure.
Les données: les données brutes (non transformées!) doivent être présentées.
Barres d’erreurs, intervalles de confiance et/ou de prédiction: pour un graphe à barres, ajouter les barres d’erreur afin de représenter la variabilité associée aux données. Pour les figures de régression, il est important d’ajouter l’intervalle de confiance ou de prédiction (dans le cas des résultats d’un modèle linéaire généralisé binomial).
Légende dans la figure: une petite légende peut être nécessaire pour distinguer les traitements (couleur, type de ligne, etc.)
Indice de significativité: on doit retrouver la valeur p d’une relation, des astérisques au-dessus de graphes à barres pour indiquer si la relation est significative, le R2 d’un modèle.
Les tableaux sont un bon choix pour présenter de l’information numérique détaillée. Ils présentent habituellement des résultats plus complexes qui seraient trop encombrants à inclure dans une figure ou dans le texte. Généralement, si les données ne peuvent être présentées en une ou deux phrases, un tableau est nécessaire. Les lignes et les colonnes doivent contenir le nom de la variable ainsi que l’unité de mesure. Un tableau résumé des statistiques peut inclure, par exemple, la moyenne, l’écart-type, les intervalles de confiances, les degrés de liberté, la valeur p et autres statistiques (comme la F value).
Légende d’une figure ou d’un tableau
Les légendes servent à compléter l’information qui est présentée dans le tableau ou la figure. On y retrouve entre autre les tailles d’échantillons (n), valeur p, des descriptions des abréviations, la méthode de collection, le nombre de réplicats, etc.
Annexe ou pas?
Les données ou les figures qui ne contribuent pas directement à l’histoire principale de votre rapport peuvent être rassemblées en annexes. Dans cette section, on peut retrouver, par exemple, des cartes du site d’échantillonnage, du matériel utilisé pour récolter les données, des tableaux avec davantage d’informations sur les modèles (pas de capture d’écran du summary!), etc.
Coder des graphiques avec ggplot
Avant de commencer, identifiez vos variables dépendantes et indépendantes. Puis, déterminez quels sont les types de ces variables avec lesquelles vous allez travailler.
Sont-elles catégoriques (et Nominales ou ordinales), ou bien numériques (et continues ou discrètes)?
C’est ce qui déterminera le choix de type de figure pour représenter vos données - un histogramme, un scatterplot, un boxplot, un barplot ? Autre?
Consultez ce lien pour mieux comprendre:
https://statisticsbyjim.com/basics/data-types/
summary(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100
## 1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300
## Median :5.800 Median :3.000 Median :4.350 Median :1.300
## Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199
## 3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800
## Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500
## Species
## setosa :50
## versicolor:50
## virginica :50
##
##
## ANOVA
La fonction de base plot() permet de visualiser rapidement un jeu de données, par exemple avec un histogramme. Par contre, vous verrez que son utilisation peut devenir limitée lorsqu’il s’agit de réaliser des figures plus complexes ou simplement modifier certains paramètres graphiques.
C’est pour cela que nous suggérons d’utiliser la fonction ggplot(), du package ggplot2. Cette fonction permet de réaliser des graphiques de façon plus intuitive et permet de les mettre en page plus facilement qu’avec la fonction plot(). Elle est aussi mieux documentée, il est donc plus facile de comprendre et utiliser les différents aspects de la fonction. Google est d’ailleurs votre meilleur allié pour la réalisation de vos figures!
Contrairement à la fonction plot(), ggplot() fonctionne par couches. Une figure ggplot() commence avec la fonction ggplot(). Elle sert à “préparer” la figure: on spécifie le jeu de données à utiliser, puis on choisit les variables qui formeront nos axes.
La fonction ggplot() nécessite deux arguments: le dataset (jeu de données), puis l’argument aes(). Ce dernier nous permet d’assigner des variables du dataset aux composantes du graphique (par exemple, les axes x et y). Voici un exemple, toujours avec le jeu de données Iris. Nous allons tester si la largeur des sépales diffère entre les espèces.
(Revoir l’atelier 3 pour l’explication de l’ANOVA, incluant les postulats et critères d’utilisation, que je passe ici)
library("ggplot2")
head(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosaboxplot(Sepal.Width ~ Species,
data = iris)
On produit d’abord le modèle linéaire, et le test d’ANOVA, avec Sepal.Width comme variable dépendante et Species comme variable indépendante.
modele.SW<-lm(Sepal.Width ~ Species, data = iris)
#Test anova:
anova(modele.SW)## Analysis of Variance Table
##
## Response: Sepal.Width
## Df Sum Sq Mean Sq F value Pr(>F)
## Species 2 11.345 5.6725 49.16 < 2.2e-16 ***
## Residuals 147 16.962 0.1154
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(modele.SW)##
## Call:
## lm(formula = Sepal.Width ~ Species, data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.128 -0.228 0.026 0.226 0.972
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.42800 0.04804 71.359 < 2e-16 ***
## Speciesversicolor -0.65800 0.06794 -9.685 < 2e-16 ***
## Speciesvirginica -0.45400 0.06794 -6.683 4.54e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3397 on 147 degrees of freedom
## Multiple R-squared: 0.4008, Adjusted R-squared: 0.3926
## F-statistic: 49.16 on 2 and 147 DF, p-value: < 2.2e-16Ensuite, on effectue le test post-hoc Tukey:
compSW <- aov(Sepal.Width ~ Species, data = iris)
TukeyHSD(compSW)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Sepal.Width ~ Species, data = iris)
##
## $Species
## diff lwr upr p adj
## versicolor-setosa -0.658 -0.81885528 -0.4971447 0.0000000
## virginica-setosa -0.454 -0.61485528 -0.2931447 0.0000000
## virginica-versicolor 0.204 0.04314472 0.3648553 0.0087802summary(compSW)## Df Sum Sq Mean Sq F value Pr(>F)
## Species 2 11.35 5.672 49.16 <2e-16 ***
## Residuals 147 16.96 0.115
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1On peut ensuite construire notre figure. D’abord, la fonction ggplot(). On spécifie le jeu de données, puis l’argument aes(), ici les axes x et y:
ggplot(iris, aes(x = Species, y = Sepal.Width))
Où sont les données? Il faut les ajouter!
C’est comme ça que ggplot() fonctionne. On crée la base de notre figure, puis on y ajoute les données à l’aide des fonctions geom_*. Par exemple, geom_points nous permet d’ajouter des données sous forme de scatterplots (donc des points), tandis que geom_line nous permet d’ajouter une ligne. Dans notre cas, nous voulons montrer nos données sous forme de boxplot, donc nous utilisons geom_boxplot. Il ne faut pas oublier d’ajouter un + après chaque fonction:
ggplot(iris, aes(x = Species, y = Sepal.Width)) +
geom_boxplot()
Pour ajouter/modifier d’autres éléments, il faut ajouter des couches. Ainsi, on peut changer les titres d’axes (labs pour labels), changer le thème de notre graphique (theme_bw est plus minimaliste), et avec la fonction theme(), retirer la grille et centrer le titre:
ggplot(iris, aes(x = Species, y = Sepal.Width)) +
geom_boxplot()+
labs(x = "Espèce",
y = "Largeur des sépales (cm)",
title = "Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5))
On peut aussi décider de changer les couleurs. Ce n’est pas toujours nécessaire ou pertinent, si ça n’ajoute pas d’information à la figure. Il faut faire attention de ne pas surcharger la figure avec le design esthétique, l’important c’est de focuser sur le message que l’on veut transmettre avec notre figure, et ne pas ajouter d’éléments qui sont distrayants. Mais pour l’exercice et se familiariser avec les paramètres graphiques, essayons:
ggplot(iris, aes(x = Species, y = Sepal.Width)) +
geom_boxplot(fill = "salmon",
color = "red")+
labs(x="Espèce",
y="Largeur des sépales (mm)",
title="Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5)) 
On a donc changé le fill (intérieur) et color (trait) de nos boxplot, à même la fonction boxplot(). On pourrait aussi changer la couleur des boîtes en fonction de l’espèce d’iris (qui est aussi notre variable x). Dans ce cas, il faut ajouter un argument aes() à la fonction boxplot(), pour assigner notre variable à un aesthetic. Ce sera notre variable Species:
#Fill:
ggplot(iris, aes(x = Species, y = Sepal.Width)) +
geom_boxplot(aes(fill = Species),
show.legend = FALSE)+
labs(x = "Espèce",
y = "Largeur des sépales (mm)",
title = "Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5)) 
#Color:
ggplot(iris, aes(x = Species, y = Sepal.Width)) +
geom_boxplot(aes(color = Species),
show.legend = FALSE)+
labs(x = "Espèce",
y = "Largeur des sépales (mm)",
title = "Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5)) 
Il y a aussi des palettes de couleurs prédéfinies dans R, qu’on peut utiliser avec la fonction scale_fill_brewer (j’ai choisi la palette “Accent”):
ggplot(iris, aes(x = Species, y = Sepal.Width)) +
geom_boxplot(aes(fill = Species),
show.legend = FALSE)+
scale_fill_brewer(palette = "Accent")+
labs(x = "Espèce",
y = "Largeur des sépales (mm)",
title = "Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5))
#Qu'est-ce qui arrive si on change Fill pour Color, dans boxplot()? Il faut aussi changer scale_fill_brewer() pour scale_color_brewer():
ggplot(iris, aes(x = Species, y=Sepal.Width)) +
geom_boxplot(aes(color = Species),
show.legend = FALSE)+
scale_color_brewer(palette = "Accent")+
labs(x = "Espèce",
y = "Largeur des sépales (mm)",
title = "Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5))
Puisque ggplot() fonctionne par couche, on peut aussi choisir d’ajouter les données brutes à notre figure. Par exemple, on ajoute avant geom_boxplot() une autre couche de données avec geom_point(), dans lequel on peut spécifier, encore une fois, l’argument aes() pour assigner la variable qui sera colorée. L’argument position définit le niveau de “jitter”, c’est-à-dire des points un peu “éparpillés” sur l’axe x. Ensuite, on pourrait ajouter l’argument alpha=0.5 dans la fonction boxplot(), pour les rendre semi-transparents:
ggplot(iris, aes(x = Species, y = Sepal.Width)) +
geom_point(aes(color = Species),
position = position_jitter(0.2),
shape = 16,
size = 2.5,
show.legend = FALSE)+
scale_color_brewer(palette = "Accent")+
geom_boxplot(alpha = 0.5)+
labs(x = "Espèce",
y = "Largeur des sépales (mm)",
title = "Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5))
On peut aussi ajouter des lettres pour montrer nos résultats de comparaison (voir https://statdoe.com/one-way-anova-and-box-plot-in-r/ pour la méthode). En gros, la fonction multcompLetters4 assigne les lettres en fonction des groupes du test de Tukey, puis on constuit une table qui contient ces lettres (LettresSW dans l’exemple), la variable Species, ainsi que la position de la lettre (en haut du quartile pour chaque boîte/valeur de Species).
Deux nouvelles libraries sont nécessaires.
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(multcompView)
compSW <- aov(Sepal.Width ~ Species, data = iris)
tukeySW <- TukeyHSD(compSW)
print(tukeySW)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Sepal.Width ~ Species, data = iris)
##
## $Species
## diff lwr upr p adj
## versicolor-setosa -0.658 -0.81885528 -0.4971447 0.0000000
## virginica-setosa -0.454 -0.61485528 -0.2931447 0.0000000
## virginica-versicolor 0.204 0.04314472 0.3648553 0.0087802cld <- multcompLetters4(compSW, tukeySW)
LettresSW <- group_by(iris, Species) %>%
summarise(mean = mean(Sepal.Width), quant = quantile(Sepal.Width, probs = 0.75)) %>%
arrange(desc(mean))
cld <- as.data.frame.list(cld$Species)
LettresSW$cld <- cld$Letters
print(LettresSW)## # A tibble: 3 × 4
## Species mean quant cld
## <fct> <dbl> <dbl> <chr>
## 1 setosa 3.43 3.68 a
## 2 virginica 2.97 3.18 b
## 3 versicolor 2.77 3 cOn ajoute nos lettres au boxplot avec la fonction geom_text(), avec comme argument aes() la table que l’on vient de créer. On peut aussi ajouter la statistique de test avec la fonction annotate(), en utilisant les coordonnées x et y dans notre figure.
ggplot(iris, aes(x = Species, y = Sepal.Width)) +
geom_point(aes(color = Species),
position = position_jitter(0.1),
shape = 16,
size = 2.5,
show.legend = FALSE)+
scale_color_brewer(palette = "Accent")+
geom_boxplot(alpha = 0.5)+
labs(x = "Espèce",
y = "Largeur des sépales (mm)",
title = "Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5))+
geom_text(data = LettresSW,
aes(x = Species, y = quant, label = cld),
size = 5,
vjust=-1, hjust =-5)+
annotate("text", x = 3, y = 4.5, label = "ANOVA, F(2)=49.16, p<0.001")
Pour bien faire, on pourrait même changer les étiquettes en x. Avec la fonction scale_x_discrete, j’ajoute “I.” pour iris à chaque étiquette et je les mets en italique, puisque ce sont des noms latins, avec l’argument axis.text.x dans la fonction theme().
ggplot(iris, aes(x = Species, y = Sepal.Width)) +
geom_point(aes(color = Species),
position = position_jitter(0.1),
shape = 16,
size = 2.5, show.legend = FALSE)+
scale_color_brewer(palette = "Accent")+
geom_boxplot(alpha = 0.5)+
labs(x = "Espèce",
y = "Largeur des sépales (mm)",
title = "Largeur des sépales chez les espèces d'iris")+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(face="italic"))+
geom_text(data = LettresSW,
aes(x = Species, y = quant, label = cld),
size = 5,
vjust=-1, hjust =-5)+
annotate("text", x = 3, y = 4.5, label = "ANOVA, F(2)=49.16, p<0.001")+
scale_x_discrete(labels = c("I. setosa", "I. versicolor", "I. virginica"))
GLM
Essayons avec un type d’analyse différent, les GLM binomiaux. Nous allons utiliser le jeu de données “survie”. Encore une fois, référez-vous à l’atelier sur les modèles GLM pour bien comprendre les analyses, j’y vais plus grossièrement ici.
Nous aurons besoin de:
library(ggeffects)D’abord, télécharger le jeu de données “survie.csv”, puis l’ouvrir dans R (*le jeu provient du cours BIO8092).
surv <- read.delim2('survie.csv', header = TRUE, sep = ';', dec = '.')
head(surv)## masse age sex survie
## 1 103.31373 60.43044 0 1
## 2 100.38447 63.11296 1 1
## 3 101.89734 60.85186 1 1
## 4 98.54894 62.10991 0 0
## 5 101.82275 59.57312 0 1
## 6 102.47745 56.40091 1 1On a donc des données de masse, d’âge et de sexe (sous forme 0/1, mâle=1 et femelle=0), ainsi que ce qui sera notre variable réponse, la survie d’une bibitte quelconque. Celle-ci est aussi notée 0/1, car évidemment, c’est une variable binaire. Puisque celle-ci ne suit pas une loi normale, nous devons utiliser un GLM.
Prenez bien le temps de comprendre comment fonctionne la régression logistique. L’interprétation des coefficients n’est pas aussi directe ou intuitive que pour la régression linéaire. Comme vous vous souvenez sûrement, les paramètres sont dans une échelle différente que l’échelle de la “réalité”: la pente nous donne plutôt le changement dans le logit de la probabilité, pour un changement d’une unité de la variable explicative.
Rappel: Les GLM fonctionnent avec ce qu’on appelle une fonction de lien (“logit”, dans le cas d’une variable binaire - il en existe d’autre pour les autres distributions, comme celle de poisson ou gamma). Cette transformatiom permet de passer de l’échelle des valeurs observées (soit 0 et 1), à une échelle qui est linéaire, qui nous permet de “linéariser” un modèle à variable réponse binaire. Ainsi: Les valeurs observées de y (0 et 1) sont transformées avec la fonction de lien f() pour être sur la même échelle que le prédicteur linéaire (a+bxi). Les paramètres du modèle, nécessairement, sont dans cette échelle transformée. C’est pour cela qu’on ne peut pas calculer directement les valeurs prédites, comme avec la régression linéaire. Plutôt, les valeurs prédites par le modèle sont obtenues en appliquant l’inverse de la transformation, donc l’inverse de la fonction de lien (f−1), au prédicteur linéaire. Ainsi, Probabilité(succès)= logit−1(a+bx).
On peut convertir les paramètres du modèle logistique en ratios des cotes, en calculant e^-b, b étant la valeur d’un paramètre donné dans un modèle (voir https://statisticsbyjim.com/probability/odds-ratio/). Ce ratio des cotes nous indique le changement de cote associé à l’augmentation d’une unité de notre variable indépendante, par exemple, l’âge (s’il s’agit d’une autre variable binaire, par exemple le sexe, l’augmentation d’unité signifie qu’on passe de zéro, soit les femelles, à 1, soit les mâles − l’ordonnée à l’origine représente le 0 de toutes les variables).
En gros,
Si le ratio des cotes est >1, la probabilité de succès (dans notre cas, la survie) augmente avec l’augmentation de notre variable dépendante.
Si le ratio des cotes est <1, la probabilité de succès (dans notre cas, la survie) diminue avec l’augmentation de notre variable dépendante.
Calculer le ratio des cotes e^-b nous permet d’avoir une idée générale de l’effet d’une variable indépendante sur notre variable réponse.
Mais ce qui vous intéresse dans le cadre de cet atelier, c’est de visualiser les effets desdites variables!
Pour simplifier les calculs, on peut utiliser la fonction ggpredict() du package ggeffects. Celle-ci nous calcule les probabilités de survie (ou de succès) ainsi qu’un intervalle de confiance. On peut ensuite utiliser ces probabilités pour réaliser une figure avec ggplot().
Reprenons notre jeu de données.
head(surv)## masse age sex survie
## 1 103.31373 60.43044 0 1
## 2 100.38447 63.11296 1 1
## 3 101.89734 60.85186 1 1
## 4 98.54894 62.10991 0 0
## 5 101.82275 59.57312 0 1
## 6 102.47745 56.40091 1 1On peut commencer avec la question suivante: quel est l’effet de la masse sur la probabilité de survie?
Pour faire notre modèle, on utilise la fonction glm, qui ressemble à la fonction lm que nous avons utilisée plus tôt. Avec glm, on ajoute l’argument family = binomial, qui spécifie la distribution avec laquelle nous travaillons, soit la distribution binomiale.
modele.masse <- with(surv, glm(survie ~ masse, family = binomial))
summary(modele.masse)##
## Call:
## glm(formula = survie ~ masse, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0758 0.2901 0.3887 0.5126 1.4173
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -45.9643 21.0574 -2.183 0.0290 *
## masse 0.4726 0.2095 2.256 0.0241 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 39.880 on 47 degrees of freedom
## Residual deviance: 34.075 on 46 degrees of freedom
## AIC: 38.075
##
## Number of Fisher Scoring iterations: 5#Quel serait le ratio des cotes?
b.masse<-modele.masse$coefficients['masse']
exp(b.masse)## masse
## 1.604189La masse aurait donc un effet positif sur la probabilité de survie. Essayons de visualiser ce résultat.
J’utilise la fonction ggpredict du package ggeffects. On spéficie d’abord le modèle à prédire, ensuite la variable que l’on souhaite visualiser en x.
predictions.masse<-ggpredict(modele.masse, "masse")## Data were 'prettified'. Consider using `terms="masse [all]"` to get
## smooth plots.predictions.masse## # Predicted probabilities of survie
##
## masse | Predicted | 95% CI
## --------------------------------
## 96 | 0.36 | [0.07, 0.81]
## 97 | 0.47 | [0.14, 0.83]
## 99 | 0.70 | [0.44, 0.87]
## 100 | 0.79 | [0.60, 0.90]
## 101 | 0.85 | [0.71, 0.93]
## 102 | 0.90 | [0.76, 0.96]
## 103 | 0.94 | [0.80, 0.98]
## 106 | 0.98 | [0.85, 1.00]##
## Not all rows are shown in the ouput. Use `print(..., n = Inf)` to show
## all rows.Dans le récapitulatif, on voit l’étendue des valeurs de masse, ainsi que la probabilité de survie prédite pour chaque valeur de masse. De plus, un intervalle de confiance à 95% a été calculé. On a donc ce qu’il faut pour réaliser notre figure.
Nous utiliserons comme source de données, la table de prédictions que nous venons de sortir.
plot.masse <- ggplot()+
geom_line(data = predictions.masse,
aes(x = x, y = predicted))
plot.masse
On voit que la valeur maximale de l’axe y est de 1, qui correspondrait à une probabilité de 100%. Il est important d’ajouter notre invervalle de confiance, avec geom_ribbon. On utilise à nouveau notre table de prédictions, puisqu’elle contient les intervalles de confiance. Avec aes(ymin = conf.low, ymax = conf.high), on délimite les limites haut et bas de notre ruban. Ces limites correspondent aux valeurs d’intervalle calculées. Avec alpha = .2, le ruban est plus translucide:
plot.masse <- ggplot()+
geom_line(data = predictions.masse,
aes(x = x, y = predicted))+
geom_ribbon(data = predictions.masse,
show.legend = FALSE,
aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high),
alpha = .2)
plot.masse
C’est bien, mais notre figure est pas mal tout nue!
On peut identifier nos axes + titre, ajouter le theme_bw(), enlever la grille + centrer le titre:
plot.masse <- ggplot()+
geom_line(data = predictions.masse,
aes(x = x, y = predicted))+
geom_ribbon(data = predictions.masse,
show.legend = FALSE,
aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high),
alpha = .2)+
labs(x = "Masse (g)",
y = "Probabilité de survie",
title = "Probabilité de survie de (?) en fonction de la masse")+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5))
plot.masse
Les valeurs de notre axe y sont de 0 à 1. On peut transformer notre axe pour qu’il soit exprimé en % sur 100, avec scale_y_continuous. De plus, ggplot() a choisi de commencer notre axe un peu en haut de zéro (même si c’est moins évident sur cette figure - retournez voir celle sans le ruban) car il optimise l’espace que prend nos données sur le graphique. Il serait plus juste de commencer notre axe à zéro. J’ajoute l’argument limits = c(0, 1) dans scale_y_continuous.
plot.masse <- ggplot()+
geom_line(data = predictions.masse,
aes(x = x, y = predicted))+
geom_ribbon(data = predictions.masse,
show.legend = FALSE,
aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high),
alpha = .2)+
labs(x = "Masse (g)",
y = "Probabilité de survie",
title = "Probabilité de survie de (?) en fonction de la masse")+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5))+
scale_y_continuous(labels = scales::percent,
limits = c(0, 1))
plot.masse
Si on voulait avoir moins de ticks sur l’axe des y (peut rendre la figure plus facilement lisible si elle est en petit format), on peut ajouter breaks = seq(0, 1, by = 0.5) dans scale_y_continuous, pour que les ticks soient à intervalle de 0.5:
plot.masse <- ggplot()+
geom_line(data = predictions.masse,
aes(x = x, y = predicted))+
geom_ribbon(data = predictions.masse,
show.legend = FALSE,
aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high),
alpha = .2)+
labs(x = "Masse (g)",
y = "Probabilité de survie",
title = "Probabilité de survie de (?) en fonction de la masse")+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5))+
scale_y_continuous(labels = scales::percent,
limits = c(0, 1),
breaks = seq(0, 1, by = 0.5))
plot.masse
Terminons la figure avec une jolie couleur, tiens: Dans geom_line(), j’ajoute color=group à l’argument aes(). On indique ainsi que le paramètre graphique de la couleur du trait sera associée aux “groupes”. Pour l’instant, on n’a qu’un seul groupe, donc ça peut sembler abstrait. Par contre, c’est ce qui nous permet, avec scale_color_manual(), de changer la couleur de notre ligne.
Même chose pour le ruban: dans geom_ribbon(), j’ajoute fill=group à l’argument aes(), pour assigner le remplissage du ruban aux “groupes”.
plot.masse <- ggplot()+
geom_line(data = predictions.masse,
aes(x = x, y = predicted, color = group),
show.legend = FALSE)+
scale_color_manual(values = c("#D28D32"))+
geom_ribbon(data = predictions.masse,
show.legend = FALSE,
aes(x = x, y = predicted, fill = group, ymin = conf.low, ymax = conf.high),
alpha = .2)+
scale_fill_manual(values = c("#D28D32"))+
labs(x = "Masse (g)",
y = "Probabilité de survie",
title = "Probabilité de survie de (?) en fonction de la masse")+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5))+
scale_y_continuous(labels = scales::percent,
limits = c(0, 1),
breaks = seq(0, 1, by = 0.5))
plot.masse
Remarquez que j’utilise color pour la ligne, et fill pour le ruban, comme nous avions abordé plus tôt avec les points dans l’ANOVA. En effet, on ne peut pas “fill” une ligne; en revanche, on peut “color” un ruban… Mais c’est sa bordure qui sera colorée:
plot.masse <- ggplot()+
geom_line(data = predictions.masse,
aes(x = x, y = predicted, color = group),
show.legend = FALSE)+
scale_color_manual(values = c("#D28D32"))+
geom_ribbon(data = predictions.masse,
show.legend = FALSE,
aes(x = x, y = predicted, color = group, ymin = conf.low, ymax = conf.high),
alpha = .2)+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5))
plot.masse
On pourrait aussi ajouter une couche de données contenant les données brutes, avec geom_point. Toutefois, on souhaite afficher les données brutes issues de surv, et non les résultats de notre modèle predictions.masse. Donc comme avec notre ANOVA, on assigne nos aes() dans geom_point(data=surv, aes(x=masse, y=survie)). J’utilise position = position_jitter(), car cela permet de mieux voir les points qui sont collés tous à la même valeur, soit 0 et 1 (aussi, détail, mais pour cette figure, je dois aussi ajouter height = 0.01 pour réduire le niveau de jitter appliqué en y, sinon le jitter s’énerve un peu trop et remplit tout l’espace - je crois que c’est lié au fait que notre variable n’a que deux valeurs…). De plus, vu que ça fait “déborder” les points un peu en dessous de 0 et un peu au dessus de 1, j’élargis les limites de l’axe y avec limits = c(-0.02, 1.02)) dans scale_y_continuous().
On termine en changeant la couleur des points et en augmentant la taille du trait avec linewidth=1.1 dans geom_line.
plot.masse <- ggplot()+
geom_line(data = predictions.masse,
aes(x = x, y = predicted, color = group),
show.legend = FALSE,
linewidth = 1.1)+
scale_color_manual(values = c("#D28D32"))+
geom_ribbon(data = predictions.masse,
show.legend = FALSE,
aes(x = x, y = predicted, fill = group, ymin = conf.low, ymax = conf.high),
alpha = .2)+
scale_fill_manual(values = c("#D28D32"))+
geom_point(data = surv,
aes(x = masse, y = survie),
alpha = .4,
color = "cornsilk4",
position = position_jitter(0.1, height = 0.01))+
labs(x = "Masse (g)",
y = "Probabilité de survie",
title = "Probabilité de survie de (?) en fonction de la masse")+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5))+
scale_y_continuous(labels = scales::percent,
limits = c(-0.02, 1.02))
plot.masse
Dans le jeu de données, on retrouve aussi la variable sexe:
head(surv)## masse age sex survie
## 1 103.31373 60.43044 0 1
## 2 100.38447 63.11296 1 1
## 3 101.89734 60.85186 1 1
## 4 98.54894 62.10991 0 0
## 5 101.82275 59.57312 0 1
## 6 102.47745 56.40091 1 1On pourrait donc faire un deuxième modèle qui inclut l’effet du sexe de la bibitte dans sa probabilité de survie.
modele.masse.sexe <- with(surv, glm(survie ~ masse + sex, family = binomial))
summary(modele.masse.sexe)##
## Call:
## glm(formula = survie ~ masse + sex, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3327 0.1828 0.3944 0.5588 1.3727
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -60.3269 26.0796 -2.313 0.0207 *
## masse 0.6094 0.2574 2.368 0.0179 *
## sex 1.3184 1.0860 1.214 0.2247
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 39.880 on 47 degrees of freedom
## Residual deviance: 32.401 on 45 degrees of freedom
## AIC: 38.401
##
## Number of Fisher Scoring iterations: 6Il ne semble pas avoir d’effet significatif du sexe sur la probabilité de survie.
Mais, pour l’exercice, faisons comme si l’effet était significatif et qu’on voulait l’illustrer dans une figure!
On calcule d’abord nos valeurs prédites de survie:
predictions.masse.sexe.test <- ggpredict(modele.masse.sexe, c("masse"))## Data were 'prettified'. Consider using `terms="masse [all]"` to get
## smooth plots.predictions.masse.sexe.test## # Predicted probabilities of survie
##
## masse | Predicted | 95% CI
## --------------------------------
## 96 | 0.24 | [0.03, 0.76]
## 97 | 0.36 | [0.08, 0.78]
## 99 | 0.66 | [0.40, 0.85]
## 100 | 0.78 | [0.59, 0.90]
## 101 | 0.87 | [0.71, 0.95]
## 102 | 0.92 | [0.77, 0.98]
## 103 | 0.96 | [0.80, 0.99]
## 106 | 0.99 | [0.87, 1.00]
##
## Adjusted for:
## * sex = 0.50##
## Not all rows are shown in the ouput. Use `print(..., n = Inf)` to show
## all rows.Remarquez dans le récapitulatif: Adjusted for: sex = 0.50
Rappelez-vous comment est construit un modèle GLM: 
Les différents termes (paramètres β)*(la valeur de x) sont additionnés pour obtenir la variable réponse.
Donc, quand on lui demande de calculer la probabilité de survie en fonction de la masse, la fonction ggpredit() utilise par défaut la valeur moyenne de l’autre variable. Il set donc la valeur de sexe à 0,5. On obtient comme figure:
plot.masse.sexe.test <- ggplot()+
geom_line(data = predictions.masse.sexe.test,
aes(x = x, y = predicted, color = group),
show.legend = FALSE,
linewidth=1.1)+
labs(title = "Probabilité de survie ajustée pour sexe = 0.5")+
theme(plot.title = element_text(hjust = 0.5))
plot.masse.sexe.test
Dans notre cas, par contre, il serait plus intéressant de pouvour visualiser cette différence.
Pour cela, il faut inclure la variable sexe dans la fonction ggpredict:
predictions.masse.sexe <- ggpredict(modele.masse.sexe, c("masse","sex"))## Data were 'prettified'. Consider using `terms="masse [all]"` to get
## smooth plots.predictions.masse.sexe## # Predicted probabilities of survie
##
## # sex = 0
##
## masse | Predicted | 95% CI
## --------------------------------
## 96 | 0.14 | [0.01, 0.75]
## 98 | 0.35 | [0.07, 0.80]
## 99 | 0.50 | [0.17, 0.84]
## 101 | 0.77 | [0.52, 0.91]
## 103 | 0.92 | [0.73, 0.98]
## 106 | 0.99 | [0.84, 1.00]
##
## # sex = 1
##
## masse | Predicted | 95% CI
## --------------------------------
## 96 | 0.38 | [0.07, 0.83]
## 98 | 0.67 | [0.31, 0.90]
## 99 | 0.79 | [0.48, 0.94]
## 101 | 0.93 | [0.70, 0.99]
## 103 | 0.98 | [0.79, 1.00]
## 106 | 1.00 | [0.86, 1.00]##
## Not all rows are shown in the ouput. Use `print(..., n = Inf)` to show
## all rows.Cette fois-ci, nous obtenons deux tables de valeurs prédites - une pour sex = 0 (femelles) et une pour sex = 1 (mâles).
Dans le précédent modèle, nous avions défini l’aes(color=group), pour colorer notre ligne. La notion de groupe était plus abstraite, car il n’y avait en réalité qu’un seul groupe. Par contre, ici, dans plot.masse.sexe, nous en avons deux - les femelles et les mâles. Donc, quand on assigne color aux group, nous colorons les deux groupes de couleurs différentes:
plot.masse.sexe <- ggplot()+
geom_line(data = predictions.masse.sexe,
aes(x = x, y = predicted, color = group),
show.legend = TRUE,
linewidth=1.1)
plot.masse.sexe
On doit ajouter nos intervalles de confiance:
plot.masse.sexe <- ggplot()+
geom_line(data = predictions.masse.sexe,
aes(x = x, y = predicted, color = group),
show.legend = TRUE)+
geom_ribbon(data = predictions.masse.sexe,
show.legend = FALSE,
aes(x = x, y = predicted, fill = group, ymin = conf.low, ymax = conf.high),
alpha = 0.2)
plot.masse.sexe
La légende doit aussi être modifiée. scale_color_manual traite la couleur du trait formé par geom_line(), et scale_fill_manual, le remplissage du ruban geom_ribbon. En même temps, avec ces mêmes fonctions, on peut spécifier le titre de la légende (on n’écrit pas “Légende”!), le nom des groupes (la variable sex est codée 0 / 1, donc on doit ajouter manuellement ces noms). On peut aussi changer les couleurs.
plot.masse.sexe <- ggplot()+
geom_line(data = predictions.masse.sexe,
aes(x = x, y = predicted, color = group),
show.legend = TRUE)+
geom_ribbon(data = predictions.masse.sexe,
show.legend = TRUE,
aes(x = x, y = predicted, fill = group, ymin = conf.low, ymax = conf.high),
alpha = 0.2)+
scale_color_manual(name = "Sexe",
labels = c("Femelle", "Mâles"),
values = c("0" = "#D28D32", "1" = "darkblue"))+
scale_fill_manual(name = "Sexe",
labels = c("Femelle", "Mâles"),
values = c("0" = "#D28D32", "1" = "darkblue"))
plot.masse.sexe
Il est important que les noms et les aes() soient constants pour chaque ligne de notre figure ggplot(). Sinon, des légendes supplémentaires sont créées:
plot.masse.sexe <- ggplot()+
geom_line(data = predictions.masse.sexe,
aes(x = x, y = predicted, color = group),
show.legend = TRUE)+
geom_ribbon(data = predictions.masse.sexe,
show.legend = TRUE,
aes(x = x, y = predicted, fill = group, ymin = conf.low, ymax = conf.high),
alpha = 0.2)+
scale_color_manual(name = "Sexe",
labels = c("Femelle", "Mâles"),
values = c("0" = "#D28D32", "1" = "darkblue"))+
scale_fill_manual(name = "Coucou",
labels = c("Femelle", "Mâles"),
values = c("0" = "#D28D32", "1" = "darkblue"))
plot.masse.sexe
On peut ajouter les éléments graphiques que nous avions codés pour la figure de notre premier modèle: geom_line(linewidth=1.1) labs() theme_bw() theme() scale_y_continuous(labels = scales::percent, limits = c(-0.02, 1.02))
plot.masse.sexe <- ggplot()+
geom_line(data = predictions.masse.sexe, aes(x=x, y = predicted, color = group),
show.legend = TRUE,
linewidth = 1.1)+
geom_ribbon(data = predictions.masse.sexe,
show.legend = TRUE,
aes(x = x, y = predicted, fill = group, ymin = conf.low, ymax = conf.high),
alpha = 0.2)+
scale_color_manual(name = "Sexe",
labels = c("Femelle", "Mâles"),
values = c("0" = "#D28D32", "1" = "darkblue"))+
scale_fill_manual(name = "Sexe",
labels = c("Femelle", "Mâles"),
values = c("0" = "#D28D32", "1" = "darkblue"))+
labs(x = "Masse (g)",
y = "Probabilité de survie",
title = "Probabilité de survie de (?) en fonction de la masse, pour les individus mâles et femelles")+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5))+
scale_y_continuous(labels = scales::percent,
limits = c(-0.02, 1.02))
plot.masse.sexe
Le titre est devenu plus long et pourrait déborder de la figure (cela dépendra des paramètres choisis au moment d’exporter la figure). On peut utiliser \n pour indiquer précisément l’endroit où couper notre texte:
plot.masse.sexe <- ggplot()+
geom_line(data = predictions.masse.sexe,
aes(x = x, y = predicted, color = group),
show.legend = TRUE,
linewidth = 1.1)+
geom_ribbon(data = predictions.masse.sexe,
show.legend = TRUE,
aes(x = x, y = predicted, fill = group, ymin = conf.low, ymax = conf.high),
alpha = 0.2)+
scale_color_manual(name = "Sexe",
labels = c("Femelle", "Mâles"),
values = c("0" = "#D28D32", "1" = "darkblue"))+
scale_fill_manual(name = "Sexe",
labels = c("Femelle", "Mâles"),
values = c("0" = "#D28D32", "1" = "darkblue"))+
labs(x = "Masse (g)",
y = "Probabilité de survie",
title = "Probabilité de survie de (?) en fonction de la\nmasse, pour les individus mâles et femelles")+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5))+
scale_y_continuous(labels = scales::percent,
limits = c(-0.02, 1.02))
plot.masse.sexe
On peut ajouter nos données brutes pour terminer la figure.
Je dois faire une dernière étape pour afficher correctement les données brutes.
surv$sex<-as.factor(surv$sex)Sinon, il traite sex comme étant une variable continue (ce qu’elle est), et on ne peut pas appliquer les paramètres graphiques correctement!
(Idéalement, faire cela avant l’analyse…!)
plot.masse.sexe <- ggplot()+
geom_line(data = predictions.masse.sexe,
aes(x = x, y = predicted, color = group),
show.legend = TRUE,
linewidth = 1.1)+
scale_color_manual(name = "Sexe",
labels = c("Femelle", "Mâles"),
values = c("0" = "#D28D32", "1" = "darkblue"))+
geom_ribbon(data = predictions.masse.sexe,
show.legend = TRUE,
aes(x = x, y = predicted, fill = group, ymin = conf.low, ymax = conf.high),
alpha = 0.2)+
scale_fill_manual(name = "Sexe",
labels = c("Femelle", "Mâles"),
values = c("0" = "#D28D32", "1" = "darkblue"))+
geom_point(data = surv, aes(x = masse, y = survie, col = sex),
alpha = .4,
position = position_jitter(0.1, height = 0.01))+
labs(x = "Masse (g)",
y = "Probabilité de survie",
title = "Probabilité de survie de (?) en fonction de la\nmasse, pour les individus mâles et femelles")+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5))+
scale_y_continuous(labels = scales::percent,
limits = c(-0.02, 1.02))+
scale_x_continuous(breaks = seq(96, 106, by = 2))
plot.masse.sexe
Si on avait une troisième variable dans notre modèle, par exemple l’âge:
head(surv)## masse age sex survie
## 1 103.31373 60.43044 0 1
## 2 100.38447 63.11296 1 1
## 3 101.89734 60.85186 1 1
## 4 98.54894 62.10991 0 0
## 5 101.82275 59.57312 0 1
## 6 102.47745 56.40091 1 1On peut l’ajouter à notre modèle:
modele.masse.sexe.age <- with(surv, glm(survie ~ masse + sex + age, family = binomial))
summary(modele.masse.sexe.age)##
## Call:
## glm(formula = survie ~ masse + sex + age, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4538 0.2002 0.3538 0.5734 1.2288
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -69.8689 30.8379 -2.266 0.0235 *
## masse 0.6183 0.2599 2.379 0.0173 *
## sex1 1.4678 1.1074 1.325 0.1850
## age 0.1434 0.2213 0.648 0.5171
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 39.880 on 47 degrees of freedom
## Residual deviance: 31.972 on 44 degrees of freedom
## AIC: 39.972
##
## Number of Fisher Scoring iterations: 6Ok, cette variable n’est pas significative non plus. Mais encore une fois, essayons de la visualiser pour l’exercice.
On débute avec le calcul des prédictions, en ajoutant la variable age:
predictions.masse.sexe.age <- ggpredict(modele.masse.sexe.age, c("masse","sex","age"))## Data were 'prettified'. Consider using `terms="masse [all]"` to get
## smooth plots.predictions.masse.sexe.age## # Predicted probabilities of survie
##
## # sex = 0
## # age = 57.64
##
## masse | Predicted | 95% CI
## --------------------------------
## 96 | 0.10 | [0.00, 0.72]
## 99 | 0.40 | [0.08, 0.83]
## 101 | 0.70 | [0.33, 0.92]
## 106 | 0.98 | [0.76, 1.00]
##
## # sex = 1
## # age = 57.64
##
## masse | Predicted | 95% CI
## --------------------------------
## 96 | 0.32 | [0.05, 0.82]
## 99 | 0.75 | [0.38, 0.93]
## 101 | 0.91 | [0.63, 0.98]
## 106 | 1.00 | [0.84, 1.00]
##
## # sex = 0
## # age = 59.98
##
## masse | Predicted | 95% CI
## --------------------------------
## 96 | 0.13 | [0.01, 0.73]
## 99 | 0.49 | [0.16, 0.83]
## 101 | 0.77 | [0.51, 0.91]
## 106 | 0.99 | [0.83, 1.00]
##
## # sex = 1
## # age = 59.98
##
## masse | Predicted | 95% CI
## --------------------------------
## 96 | 0.39 | [0.07, 0.84]
## 99 | 0.80 | [0.49, 0.95]
## 101 | 0.93 | [0.71, 0.99]
## 106 | 1.00 | [0.87, 1.00]
##
## # sex = 0
## # age = 62.32
##
## masse | Predicted | 95% CI
## --------------------------------
## 96 | 0.17 | [0.01, 0.80]
## 99 | 0.57 | [0.18, 0.89]
## 101 | 0.82 | [0.51, 0.95]
## 106 | 0.99 | [0.84, 1.00]
##
## # sex = 1
## # age = 62.32
##
## masse | Predicted | 95% CI
## --------------------------------
## 96 | 0.47 | [0.07, 0.91]
## 99 | 0.85 | [0.45, 0.98]
## 101 | 0.95 | [0.68, 0.99]
## 106 | 1.00 | [0.87, 1.00]##
## Not all rows are shown in the ouput. Use `print(..., n = Inf)` to show
## all rows.Nous obtenons 6 tables de prédictions. Puisque age est une variable continue, ggpredict nous a choisi 3 valeurs de age représentatives de la distribution et a calculé un y prédit pour chacunes de ces valeurs et pour les 2 niveaux de sex.
Ici nous avons deux options. Puisque nous avons une troisième variable indépendante à représenter, on peut lui assigner un aes(), par exemple linetype.
Nous avons donc, dans nos tables de prédiction:
masse = x = x
survie = y = predicted
sex = color = group
Dans la table de prédiction, bien que ce ne soit pas affiché, la troisième variable indépendante est appelée facet, comme on peut remarquer dans
str(predictions.masse.sexe.age)## Classes 'ggeffects' and 'data.frame': 66 obs. of 7 variables:
## $ x : int 96 96 96 96 96 96 97 97 97 97 ...
## $ predicted: num 0.0958 0.1291 0.1717 0.315 0.3914 ...
## $ std.error: num 1.63 1.48 1.51 1.16 1.08 ...
## $ conf.low : num 0.00433 0.00803 0.01063 0.04535 0.07209 ...
## $ conf.high: num 0.721 0.731 0.8 0.817 0.842 ...
## $ group : Factor w/ 2 levels "0","1": 1 1 1 2 2 2 1 1 1 2 ...
## $ facet : Factor w/ 3 levels "57.64","59.98",..: 1 2 3 1 2 3 1 2 3 1 ...
## - attr(*, "numeric.facet")= logi TRUE
## - attr(*, "legend.labels")= chr [1:2] "0" "1"
## - attr(*, "x.is.factor")= chr "0"
## - attr(*, "continuous.group")= logi FALSE
## - attr(*, "rawdata")='data.frame': 48 obs. of 5 variables:
## ..$ response: int [1:48] 1 1 1 0 1 1 1 1 1 1 ...
## ..$ x : num [1:48] 103.3 100.4 101.9 98.5 101.8 ...
## ..$ group : Factor w/ 2 levels "0","1": 1 2 2 1 1 2 1 2 1 2 ...
## ..$ facet : num [1:48] 60.4 63.1 60.9 62.1 59.6 ...
## ..$ rowname : chr [1:48] "1" "2" "3" "4" ...
## - attr(*, "title")= chr "Predicted probabilities of survie"
## - attr(*, "x.title")= chr "masse"
## - attr(*, "y.title")= chr "survie"
## - attr(*, "legend.title")= chr "sex"
## - attr(*, "constant.values")= Named list()
## - attr(*, "terms")= chr [1:3] "masse" "sex" "age"
## - attr(*, "original.terms")= chr [1:3] "masse" "sex" "age"
## - attr(*, "at.list")=List of 3
## ..$ masse: int [1:11] 96 97 98 99 100 101 102 103 104 105 ...
## ..$ sex : chr [1:2] "0" "1"
## ..$ age : num [1:3] 57.6 60 62.3
## - attr(*, "ci.lvl")= num 0.95
## - attr(*, "type")= chr "fe"
## - attr(*, "response.name")= chr "survie"
## - attr(*, "back.transform")= logi TRUE
## - attr(*, "response.transform")= chr "survie"
## - attr(*, "untransformed.predictions")= num [1:66] 0.0958 0.1291 0.1717 0.315 0.3914 ...
## - attr(*, "family")= chr "binomial"
## - attr(*, "link")= chr "logit"
## - attr(*, "logistic")= chr "1"
## - attr(*, "link_inverse")=function (eta)
## - attr(*, "link_function")=function (mu)
## - attr(*, "is.trial")= chr "0"
## - attr(*, "fitfun")= chr "glm"
## - attr(*, "verbose")= logi TRUE
## - attr(*, "model.name")= chr "modele.masse.sexe.age"Dans ce cas, on ajoute l’aes() linetype à geom_line assigné à facet:
plot.masse.sexe.age1 <- ggplot()+
geom_line(data = predictions.masse.sexe.age,
aes(x = x, y = predicted, color = group, linetype = facet),
show.legend = TRUE,
linewidth=1.1)
plot.masse.sexe.age1
On peut voir dans la légende les 3 valeurs de age qui ont été utilisées pour le calcul.
La deuxième option serait d’utiliser des facettes. Ceci crée des graphiques côte à côte avec différents niveaux d’une variable choisie. Dans ce cas, j’enlève l’aes() qu’on vient d’ajouter, et j’ajoute la fonction facet_wrap(~facet), qui crée 3 facettes avec notre “variable” facet:
plot.masse.sexe.age2 <- ggplot()+
geom_line(data = predictions.masse.sexe.age,
aes(x = x, y = predicted, color = group),
show.legend = TRUE,
linewidth = 1.1)+
facet_wrap(~ facet)
plot.masse.sexe.age2
On peut terminer ce graphique, avec le même code que les précédents. J’ajoute labeller dans facet_wrap pour changer le titre des étiquettes. J’ai arrondi l’âge pour les étiquettes.
plot.masse.sexe.age2 <- ggplot()+
geom_line(data = predictions.masse.sexe.age,
aes(x = x, y = predicted, color = group),
show.legend = TRUE,
linewidth = 1.1)+
scale_color_manual(name = "Sexe",
labels = c("Femelle", "Mâles"),
values = c("0" = "#D28D32", "1" = "darkblue"))+
geom_ribbon(data = predictions.masse.sexe,
show.legend = TRUE,
aes(x = x, y = predicted, fill = group, ymin = conf.low, ymax = conf.high),
alpha = 0.2)+
scale_fill_manual(name = "Sexe",
labels = c("Femelle", "Mâles"),
values = c("0" = "#D28D32", "1" = "darkblue"))+
geom_point(data = surv,
aes(x = masse, y = survie, col = sex),
alpha = .4,
position = position_jitter(0.1, height = 0.01))+
labs(x = "Masse (g)",
y = "Probabilité de survie",
title = "Probabilité de survie de (?) en fonction de la\nmasse et de l'âge, pour les individus mâles et femelles")+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5))+
scale_y_continuous(labels = scales::percent,
limits = c(-0.02, 1.02))+
scale_x_continuous(breaks = seq(96, 106, by = 2))+
facet_wrap(~ facet, labeller = labeller(facet = c( "57.64" = "58 jours", "59.98" = "60 jours", "62.32" = "62 jours")))
plot.masse.sexe.age2
Vous pouvez vous pratiquer en complétant la première option!
plot.masse.sexe.age1
Régression linéaire simple
Maintenant qu’on a fait un exemple plus compliqué, essayons d’illustrer un modèle linéaire simple.
Reprenons les données iris()
head(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosaOn peut créer un modèle linéaire pour prédire Sepal.Length en fonction de Petal.Length:
Encore une fois, ce n’est qu’un exemple pour s’exercer avec ggplot(). Il ne faut pas oublier les étapes préliminaires, que je saute ici, avant de se lancer dans un modèle.
regr.SL <- lm(Sepal.Length ~ Petal.Length, data = iris)
summary(regr.SL)##
## Call:
## lm(formula = Sepal.Length ~ Petal.Length, data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.24675 -0.29657 -0.01515 0.27676 1.00269
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.30660 0.07839 54.94 <2e-16 ***
## Petal.Length 0.40892 0.01889 21.65 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4071 on 148 degrees of freedom
## Multiple R-squared: 0.76, Adjusted R-squared: 0.7583
## F-statistic: 468.6 on 1 and 148 DF, p-value: < 2.2e-16Enfin quelque chose de significatif! Super, nous pouvons représenter cettre relation.
Avec une régression linéaire simple, évidemment, nous n’avons pas besoin d’une fonction comme ggpredict pour calculer y.
On peut commencer par spécifier les aes() du jeu de données brutes directement dans ggplot(); ensuite, nous visualisons ces données avec geom_point; puis, geom_smooth(method = lm) (calcule, puis) ajoute la courbe de la régression.
plot.regr.SL <- ggplot(data = iris, aes(x = Petal.Length, y = Sepal.Length))+
geom_point()+
geom_smooth(method = lm,
se = TRUE)
plot.regr.SL## `geom_smooth()` using formula = 'y ~ x'
On peut embellir notre graphique:
plot.regr.SL <- ggplot(data = iris, aes(x = Petal.Length, y = Sepal.Length))+
geom_point(col = "#8F1E5F")+
geom_smooth(method = lm,
se = TRUE,
col = "#2B091D",
fill = "#591130")+
labs(x = "Longueur du pétale (cm)",
y = "Longueur du sépale (mm)",
title = "Longueur du pétale en fonction de la longueur du sépale chez les iris")+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5))
plot.regr.SL## `geom_smooth()` using formula = 'y ~ x'
On retourne voir le résumé de notre modèle pour ajouter la statistique:
summary(regr.SL)##
## Call:
## lm(formula = Sepal.Length ~ Petal.Length, data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.24675 -0.29657 -0.01515 0.27676 1.00269
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.30660 0.07839 54.94 <2e-16 ***
## Petal.Length 0.40892 0.01889 21.65 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4071 on 148 degrees of freedom
## Multiple R-squared: 0.76, Adjusted R-squared: 0.7583
## F-statistic: 468.6 on 1 and 148 DF, p-value: < 2.2e-16Et on l’ajoute au graphique:
plot.regr.SL <- ggplot(data = iris, aes(x = Petal.Length, y = Sepal.Length))+
geom_point(col="#8F1E5F")+
geom_smooth(method = lm,
se = TRUE, col = "#2B091D", fill = "#591130")+
labs(x = "Longueur du pétale (cm)",
y = "Longueur du sépale (mm)",
title = "Longueur du pétale en fonction de la longueur du sépale chez les iris")+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5))+
annotate("text",
x = 2.5, y = 7.75,
label = "R2 = 0.76, F(1, 148) = 468.6, p = < 0.001")
plot.regr.SL## `geom_smooth()` using formula = 'y ~ x'
C’est important d’ajouter un intervalle de confiance, ainsi qu’un intervalle de prédiction. L’intervalle de confiance s’affiche avec geom_smooth(se = TRUE). Pour ce qui est de l’intervalle de prédiction, la fonction predict() calcule ces valeurs, puis colle ce nouveau dataframe à iris pour créer iris.int.pred. Nous utiliserons maintenant cette nouvelle table pour notre régression.
int.pred <- predict(regr.SL, interval = "prediction")## Warning in predict.lm(regr.SL, interval = "prediction"): predictions on current data refer to _future_ responsesiris.int.pred <- cbind(iris, int.pred)
head(iris.int.pred)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species fit lwr
## 1 5.1 3.5 1.4 0.2 setosa 4.879095 4.067202
## 2 4.9 3.0 1.4 0.2 setosa 4.879095 4.067202
## 3 4.7 3.2 1.3 0.2 setosa 4.838202 4.025897
## 4 4.6 3.1 1.5 0.2 setosa 4.919987 4.108491
## 5 5.0 3.6 1.4 0.2 setosa 4.879095 4.067202
## 6 5.4 3.9 1.7 0.4 setosa 5.001771 4.191017
## upr
## 1 5.690987
## 2 5.690987
## 3 5.650508
## 4 5.731483
## 5 5.690987
## 6 5.812526Dans geom_ribbon(aes(ymin = lwr, ymax = upr), les aes() permettent d’assigner les bornes supérieures et inférieures (en y) du ruban aux variables lwr et upr que nous venons de calculer et ajouter à notre table.
plot.regr.SL <- ggplot(data = iris.int.pred, aes(x = Petal.Length, y = Sepal.Length))+
geom_point(col = "#8F1E5F")+
geom_ribbon(aes(ymin = lwr, ymax = upr),
alpha = .2,
fill = "#591130")+
geom_smooth(method = lm,
se = TRUE,
col = "#2B091D",
fill = "#591130")+
labs(x = "Longueur du pétale (cm)",
y = "Longueur du sépale (mm)",
title = "Longueur du pétale en fonction de la longueur du sépale chez les iris")+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(hjust = 0.5))+
annotate("text",
x = 2.5, y = 7.75,
label = "R2 = 0.76, F(1, 148) = 468.6, p = < 0.001")
plot.regr.SL## `geom_smooth()` using formula = 'y ~ x'
Wow!
Cheat sheet ggplot()
Maintenant que vous êtes plus familiers avec ggplot(), voici un lien pour une cheat sheet pour ggplot():
https://rstudio.github.io/cheatsheets/html/data-visualization.html
Elle contient toute l’information sur les différents aes(), geom_*(), scale_*, ainsi que sur les thèmes, les ajustements de position, et les légendes, par exemple.
À considérer lors de la création de figures
Comme mentionné plus tôt, il faut être judicieux dans le choix de notre symbologie. En effet, les couleurs et les symboles ont souvent un sens qui nous est intuitif, et représentent ainsi de l’information - quelles couleurs choisiriez-vous pour illustrer une température chaude? Froide?
Ces choix de couleur contribuent à donner un sens à vos figures. Il existe des palettes de couleurs conçue pour illustrer des données, par exemple, d’élévation, ou de température. Certaines sont conçues pour maximiser le contraste entre les extrêmes des valeurs d’une variable. 
Il existe même une palette, “Okabe-Ito”, contenant des couleurs qui sont plus accessibles pour les gens ayant une variété de problèmes visuels. D’autres palettes sont simplement esthétiques.
On doit aussi penser au type de variable que nous souhaitons illustrer. Par exemple, il serait moins approprié d’utiliser la valeur (donc de pâle à foncé) d’une même couleur pour colorer des variables catégoriques qui ne sont pas ordonnées (pensez aux boxplots des espèces d’iris), puisque celle-ci traduit la notion d’ordre et de progression. Une palette avec des couleurs aléatoires comme celle que nous avons utilisée est mieux dans ce cas. La valeur serait toutefois appropriée pour des données de concentration, par exemple.
Toutefois, certaines publications scientifiques n’acceptent que les figures en noir et blanc; dans ce cas, on peut utiliser des symboles, des textures (par exemple, le grain peut être utilisé au lieu de la valeur), ou différentes lignes pointillées, par exemple.
De plus, assurez-vous que vos axes sont représentatifs de vos données.
Pour terminer, utilisez une fonction pour exporter vos figures. “Copier-coller” ne suffit pas, car la résolution ne sera pas assez bonne.
Voici un exemple:
ggsave(path = "votre_chemin_comme_dans_setwd()", "titre_de_votre_figures.jpeg", width = 6, height = 6, dpi=700)
Ordinations
Comme pour les analyses univariées, cet atelier est seulement axé sur la création de graphiques. Pour une introduction aux analyses multivariées (ordinations, PERMANOVA, Procrustes), consulter l’atelier 6.
Importer les données multivariées
library(vegan)## Loading required package: permute## Loading required package: lattice## This is vegan 2.6-4data(dune) # matrice de communautés
data(dune.env) # données environnementalesCe jeu de données contient une matrice d’abondance (dune) d’espèces végétales avec leur classe de couverture pour 20 sites. Le dataframe (dune.env) contient 5 variables denvironnementales.
PCA
Préparation des données
Transformation Hellinger
dune.hel <- decostand(dune, method="hellinger")Faire une PCA (avec la fonction prcomp cette fois-ci, car l’objet qu’on va créer peut être “lu” plus tard par la fonction ggord)
pca.dune <- prcomp(dune.hel)
summary(pca.dune)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 0.4122 0.3333 0.2374 0.2097 0.20242 0.17376 0.14506
## Proportion of Variance 0.3057 0.1998 0.1013 0.0791 0.07371 0.05431 0.03785
## Cumulative Proportion 0.3057 0.5054 0.6068 0.6859 0.75960 0.81391 0.85177
## PC8 PC9 PC10 PC11 PC12 PC13 PC14
## Standard deviation 0.13392 0.11745 0.10441 0.09995 0.09550 0.07985 0.07188
## Proportion of Variance 0.03226 0.02481 0.01961 0.01797 0.01641 0.01147 0.00929
## Cumulative Proportion 0.88403 0.90885 0.92846 0.94643 0.96283 0.97430 0.98360
## PC15 PC16 PC17 PC18 PC19 PC20
## Standard deviation 0.05424 0.05240 0.04063 0.03373 0.02536 1.816e-17
## Proportion of Variance 0.00529 0.00494 0.00297 0.00205 0.00116 0.000e+00
## Cumulative Proportion 0.98889 0.99383 0.99680 0.99884 1.00000 1.000e+00Dans l’atelier 6, nous avons fait nos PCA avec la fonction rda. Lorsqu’on compare le sommaire des deux analyses, on voit que les deux fonctions donnent les mêmes résultats.
pca2.dune <- rda(dune.hel)
summary(pca2.dune)##
## Call:
## rda(X = dune.hel)
##
## Partitioning of variance:
## Inertia Proportion
## Total 0.5559 1
## Unconstrained 0.5559 1
##
## Eigenvalues, and their contribution to the variance
##
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Eigenvalue 0.1699 0.1111 0.05634 0.04397 0.04097 0.03019 0.02104
## Proportion Explained 0.3057 0.1998 0.10135 0.07910 0.07371 0.05431 0.03785
## Cumulative Proportion 0.3057 0.5054 0.60679 0.68590 0.75960 0.81391 0.85177
## PC8 PC9 PC10 PC11 PC12 PC13
## Eigenvalue 0.01794 0.01379 0.01090 0.009991 0.00912 0.006376
## Proportion Explained 0.03226 0.02481 0.01961 0.017972 0.01641 0.011469
## Cumulative Proportion 0.88403 0.90885 0.92846 0.946428 0.96283 0.974303
## PC14 PC15 PC16 PC17 PC18 PC19
## Eigenvalue 0.005166 0.002942 0.002745 0.001651 0.001138 0.0006432
## Proportion Explained 0.009293 0.005292 0.004939 0.002969 0.002047 0.0011570
## Cumulative Proportion 0.983597 0.988888 0.993827 0.996796 0.998843 1.0000000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores: 1.80276
##
##
## Species scores
##
## PC1 PC2 PC3 PC4 PC5 PC6
## Achimill -0.2241355 0.057048 -0.068674 -0.173281 -0.0498964 -0.0006022
## Agrostol 0.4142446 -0.180068 0.020078 -0.020756 0.0029054 0.0193453
## Airaprae -0.0350694 0.154644 0.116317 -0.039999 -0.1384374 -0.0376930
## Alopgeni 0.1386550 -0.318444 0.186419 0.001329 0.0030171 -0.0215368
## Anthodor -0.1962910 0.238579 0.095768 -0.179434 -0.0380470 -0.0471885
## Bellpere -0.1148332 -0.054380 -0.057640 0.064957 0.0010453 0.0896219
## Bromhord -0.1410948 -0.053134 -0.049586 -0.057365 0.0210053 0.0938974
## Chenalbu 0.0121539 -0.025919 0.036102 -0.015353 -0.0038179 0.0173353
## Cirsarve -0.0026599 -0.031975 0.004872 0.021813 -0.0182940 0.0064877
## Comapalu 0.1100812 0.048635 -0.061981 -0.028257 0.0102828 0.1064664
## Eleopalu 0.3796078 0.069099 -0.185862 -0.064110 0.0144632 0.0035959
## Elymrepe -0.1456860 -0.248314 -0.094249 0.010195 -0.1482516 -0.0563573
## Empenigr 0.0004895 0.057213 0.063829 0.026314 -0.0393467 -0.0066958
## Hyporadi -0.0559119 0.197829 0.139068 0.036814 -0.1427093 -0.0326104
## Juncarti 0.2478729 -0.012497 -0.092252 0.005525 0.0265396 -0.1569375
## Juncbufo 0.0190348 -0.121457 0.179346 -0.051491 0.0931930 -0.0211675
## Lolipere -0.3366259 -0.173707 -0.206664 0.154192 0.0092155 -0.0325721
## Planlanc -0.2547353 0.197554 -0.034482 -0.039980 0.1469090 -0.0256698
## Poaprat -0.2808972 -0.142814 -0.082463 0.084918 -0.0502270 -0.0278725
## Poatriv -0.1317430 -0.363370 0.059442 -0.137436 0.0581366 -0.0136469
## Ranuflam 0.2783344 0.024399 -0.081013 -0.052662 0.0009058 0.0274683
## Rumeacet -0.1079283 -0.024094 0.054027 -0.122507 0.2108960 -0.0801798
## Sagiproc 0.0409392 -0.082141 0.246528 0.127747 -0.0148906 -0.0086761
## Salirepe 0.0651846 0.156446 0.007210 0.136628 -0.0166452 -0.0581183
## Scorautu -0.0562239 0.158686 0.094988 0.095244 0.0265794 0.1029567
## Trifprat -0.1001443 0.028383 -0.024352 -0.095019 0.1389927 -0.0417313
## Trifrepe -0.0471194 -0.006999 0.052812 0.045812 0.1395303 0.2537096
## Vicilath -0.0541878 0.052774 -0.023104 0.114929 0.0381791 0.0343033
## Bracruta 0.0852752 0.107955 -0.003783 0.181534 0.2198174 -0.1398609
## Callcusp 0.2045115 0.057120 -0.104072 -0.065216 -0.0146101 0.0579472
##
##
## Site scores (weighted sums of species scores)
##
## PC1 PC2 PC3 PC4 PC5 PC6
## 1 -0.423423 -0.42754 -6.533e-01 0.01583 -0.761705 -0.463469
## 2 -0.356333 -0.36816 -1.991e-01 -0.08370 -0.355682 0.630146
## 3 -0.080900 -0.55215 -5.753e-02 0.27275 -0.217514 0.016724
## 4 -0.041005 -0.49292 7.510e-02 0.33626 -0.282018 0.100013
## 5 -0.458003 0.05220 -1.432e-01 -0.54110 0.351891 -0.117334
## 6 -0.389342 0.20337 -1.022e-01 -0.30549 0.758789 -0.245696
## 7 -0.451813 0.06864 -6.835e-02 -0.41821 0.585536 -0.138732
## 8 0.295113 -0.29375 -6.362e-02 0.19828 0.010903 -0.185034
## 9 0.035973 -0.49264 2.714e-01 0.07945 0.086166 -0.419257
## 10 -0.453094 0.13639 -2.391e-01 -0.12201 0.144355 0.462946
## 11 -0.273214 0.29633 3.931e-05 0.87663 0.126762 0.097096
## 12 0.246622 -0.33173 9.205e-01 -0.03527 0.493526 -0.017653
## 13 0.226908 -0.48390 6.740e-01 -0.28663 -0.071279 0.323642
## 14 0.570818 0.25288 -3.165e-01 -0.32869 -0.085439 1.196310
## 15 0.654415 0.28845 -3.733e-01 0.01035 0.196968 0.002255
## 16 0.720030 -0.10689 -2.727e-01 -0.35992 0.007984 -0.504518
## 17 -0.317465 0.75272 3.395e-01 -0.64284 -0.803236 -0.262492
## 18 -0.201131 0.39818 -2.007e-01 0.89875 0.365682 0.086314
## 19 0.006263 0.73205 8.167e-01 0.33669 -0.503443 -0.085674
## 20 0.689579 0.36848 -4.077e-01 0.09888 -0.048247 -0.475587PERMANOVA. On va tester si la composition végétale diffère entre différents type d’aménagement (facteur; variable quantitative) des sites.
dune.dist <- vegdist(decostand(dune, method='hellinger'), method='euclid')
disper.dune <- betadisper(dune.dist, dune.env$Management)
anova(disper.dune)## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 3 0.17812 0.059373 1.5729 0.2349
## Residuals 16 0.60397 0.037748permanova.dune <- adonis2(dune.dist ~ Management, data = dune.env)
permanova.dune## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dune.dist ~ Management, data = dune.env)
## Df SumOfSqs R2 F Pr(>F)
## Management 3 3.1672 0.29986 2.2842 0.004 **
## Residual 16 7.3950 0.70014
## Total 19 10.5621 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Le test montre que le type d’aménagement a un effet statistiquement significatif sur la composition en espèces végétales des sites.
ggord + ggplot
On va créer un graphique pour présenter la PCA avec les différents types d’aménagement.
Extraire les scores des sites et des espèces, à partir de l’objet créer par prcomp. Cela va nous permettre de mettre éventuellement seulement les espèces désirées sur le graphique.
site.scores<-pca.dune$x
(site.scores.df <- data.frame(site.scores))## PC1 PC2 PC3 PC4 PC5
## 1 -0.42202236 -0.34450749 -3.749271e-01 0.008024023 -0.372799776
## 2 -0.35515437 -0.29665390 -1.142664e-01 -0.042440109 -0.174080679
## 3 -0.08063254 -0.44491000 -3.301696e-02 0.138292966 -0.106457407
## 4 -0.04086906 -0.39718712 4.310306e-02 0.170494078 -0.138027375
## 5 -0.45648854 0.04206103 -8.218385e-02 -0.274349726 0.172225272
## 6 -0.38805443 0.16386896 -5.863489e-02 -0.154889885 0.371372259
## 7 -0.45031849 0.05530804 -3.922469e-02 -0.212042295 0.286577531
## 8 0.29413758 -0.23669691 -3.651072e-02 0.100532017 0.005336318
## 9 0.03585360 -0.39695990 1.557600e-01 0.040285686 0.042171903
## 10 -0.45159580 0.10989979 -1.372343e-01 -0.061864199 0.070651366
## 11 -0.27231082 0.23877731 2.256124e-05 0.444472361 0.062040814
## 12 0.24580687 -0.26730006 5.282509e-01 -0.017884642 0.241545549
## 13 0.22615753 -0.38991817 3.868134e-01 -0.145331015 -0.034885901
## 14 0.56893064 0.20376783 -1.816319e-01 -0.166655594 -0.041816081
## 15 0.65225130 0.23242747 -2.142232e-01 0.005247591 0.096401575
## 16 0.71764913 -0.08613312 -1.564761e-01 -0.182487145 0.003907619
## 17 -0.31641541 0.60652639 1.948300e-01 -0.325939071 -0.393125861
## 18 -0.20046585 0.32084547 -1.151960e-01 0.455690043 0.178975130
## 19 0.00624186 0.58987086 4.687034e-01 0.170712594 -0.246398835
## 20 0.68729916 0.29691351 -2.339574e-01 0.050132321 -0.023613421
## PC6 PC7 PC8 PC9 PC10
## 1 -0.1947189103 0.14297479 0.194301690 0.112376805 0.032457596
## 2 0.2647456146 -0.17427545 -0.044272021 -0.016720326 -0.143103167
## 3 0.0070264595 -0.07648163 -0.075063932 0.010914507 0.127904407
## 4 0.0420189647 -0.17446434 -0.013889808 -0.162887863 0.140023693
## 5 -0.0492961680 -0.23585887 0.081909418 -0.045412441 0.024452056
## 6 -0.1032250470 0.07513035 0.096314316 -0.011148451 0.095001576
## 7 -0.0582862075 0.08584994 0.002778682 0.046034551 -0.050845625
## 8 -0.0777390258 0.19094461 -0.154446492 -0.153598933 -0.114991844
## 9 -0.1761440885 0.01711273 0.228755610 -0.015221320 -0.171424775
## 10 0.1944992350 -0.05049491 -0.210307980 -0.082508198 -0.058840139
## 11 0.0407933227 0.29970308 0.011721456 -0.056702121 0.077192077
## 12 -0.0074164915 -0.02536257 0.038872424 0.029358985 0.088641851
## 13 0.1359729727 0.14248707 -0.125859881 0.202329767 -0.093788726
## 14 0.5026102903 0.10186304 0.209182998 0.066090412 0.102879773
## 15 0.0009474427 0.02473974 0.061662868 -0.227259133 -0.140607818
## 16 -0.2119652616 -0.02270390 -0.177098525 -0.022669604 0.177063594
## 17 -0.1102816574 0.12212982 -0.151255687 -0.008799111 0.014706687
## 18 0.0362632716 -0.11150681 -0.122134229 0.204232852 -0.018315026
## 19 -0.0359944048 -0.15838135 0.156730498 -0.067611130 -0.008421206
## 20 -0.1998103115 -0.17340534 -0.007901405 0.199200753 -0.079984983
## PC11 PC12 PC13 PC14 PC15
## 1 -0.054588101 0.034462703 -0.035701587 0.117878451 -0.002213425
## 2 0.019457009 0.121961523 -0.072668906 -0.027268598 -0.114056583
## 3 0.189426459 0.098304109 0.102593714 0.003582509 0.063190553
## 4 -0.139198735 -0.233732456 0.058763244 -0.004413976 0.004692423
## 5 0.048675368 -0.054623058 0.056167712 -0.061130295 -0.067831090
## 6 -0.016507628 0.106475665 0.122588589 0.013653305 0.009980931
## 7 -0.118264835 -0.044205396 -0.022019777 0.012121836 0.002703392
## 8 -0.062119291 0.118226841 0.103731357 -0.096723711 -0.011208506
## 9 0.113994955 -0.061220658 -0.078065568 -0.121604909 0.085266972
## 10 -0.096890649 0.019846981 -0.109878765 0.042312047 0.134363570
## 11 -0.028975224 -0.049496469 -0.072451734 -0.063347298 -0.070411540
## 12 0.042099025 0.020074642 -0.142911945 0.072984250 -0.032959591
## 13 -0.061978791 -0.062186344 0.112696440 0.064723117 -0.016580691
## 14 0.015038111 -0.002557763 0.007135486 -0.066269669 0.032849963
## 15 0.119316583 -0.069303367 0.026613528 0.170846882 -0.023002423
## 16 -0.009598986 0.085246184 -0.110541292 -0.022926853 -0.014930180
## 17 0.148298447 -0.106347535 -0.005376556 -0.035117493 -0.003259164
## 18 0.130118499 -0.034162850 0.015833905 0.027468954 0.003240154
## 19 -0.116765892 0.158724721 0.036465495 0.028445289 0.022972462
## 20 -0.121536323 -0.045487474 0.007026659 -0.055213837 -0.002807226
## PC16 PC17 PC18 PC19 PC20
## 1 0.0127472182 -0.029903023 -0.024856595 -0.012336708 4.336809e-17
## 2 -0.0446611349 0.005305086 0.022078375 0.031007751 8.673617e-18
## 3 -0.0271552448 0.058753839 0.046914440 -0.030425313 7.979728e-17
## 4 -0.0437358342 -0.010564791 -0.011527629 0.023511185 -4.683753e-17
## 5 0.1097759923 -0.013436149 -0.016784278 -0.042864301 -2.428613e-17
## 6 -0.0306847435 -0.055202850 0.037060951 0.049957778 -1.040834e-17
## 7 -0.0816920274 0.113904666 -0.027202323 -0.012878277 -3.469447e-17
## 8 -0.0336516196 -0.042583349 -0.052191470 -0.028137562 0.000000e+00
## 9 0.0244642383 0.006081289 -0.000449554 0.027768859 7.632783e-17
## 10 0.0459636131 -0.036301400 0.020925899 -0.012264575 -5.724587e-17
## 11 0.0544619654 0.019844250 0.051357700 -0.009062018 -7.806256e-17
## 12 -0.0627871904 -0.055322352 -0.003578525 -0.032477460 5.377643e-17
## 13 0.0862150522 0.009678460 0.014880760 0.014241811 5.724587e-17
## 14 -0.0147797128 -0.003349925 -0.024719508 -0.005808041 2.775558e-17
## 15 0.0061529347 0.016725353 0.012680768 0.001031098 7.285839e-17
## 16 0.0634859914 0.037481480 -0.023923219 0.034616275 4.163336e-17
## 17 -0.0497449145 -0.020609099 0.001575928 0.004009650 8.673617e-18
## 18 -0.0006542089 -0.012827698 -0.067164565 0.018814567 -1.387779e-17
## 19 0.0295200565 0.034403175 -0.015520396 0.002854336 7.459311e-17
## 20 -0.0432404313 -0.022076960 0.060443242 -0.021559055 8.326673e-17sp.scores<-pca.dune$rotation
(sp.scores.df <- data.frame(sp.scores))## PC1 PC2 PC3 PC4 PC5
## Achimill -0.2248791555 0.07079821 -0.119661537 -0.341760328 -0.101948354
## Agrostol 0.4156189533 -0.22347029 0.034985202 -0.040936713 0.005936338
## Airaprae -0.0351857956 0.19191755 0.202676001 -0.078888514 -0.282855712
## Alopgeni 0.1391150350 -0.39519834 0.324826566 0.002621325 0.006164517
## Anthodor -0.1969422322 0.29608445 0.166870274 -0.353895610 -0.077737618
## Bellpere -0.1152141976 -0.06748736 -0.100435394 0.128114504 0.002135723
## Bromhord -0.1415629580 -0.06594128 -0.086400941 -0.113139773 0.042917985
## Chenalbu 0.0121942554 -0.03216665 0.062905676 -0.030280595 -0.007800813
## Cirsarve -0.0026687294 -0.03968196 0.008489110 0.043021061 -0.037378414
## Comapalu 0.1104463921 0.06035712 -0.107998204 -0.055730602 0.021009915
## Eleopalu 0.3808672608 0.08575416 -0.323855812 -0.126443878 0.029551309
## Elymrepe -0.1461693625 -0.30816535 -0.164224860 0.020107591 -0.302907981
## Empenigr 0.0004910767 0.07100360 0.111218624 0.051899454 -0.080393271
## Hyporadi -0.0560973565 0.24551212 0.242318916 0.072606975 -0.291583945
## Juncarti 0.2486953063 -0.01550957 -0.160744874 0.010897632 0.054225713
## Juncbufo 0.0190979141 -0.15073232 0.312500987 -0.101554267 0.190412207
## Lolipere -0.3377427168 -0.21557574 -0.360101987 0.304110194 0.018829208
## Planlanc -0.2555804426 0.24517118 -0.060083932 -0.078852590 0.300164812
## Poaprat -0.2818291593 -0.17723703 -0.143688346 0.167481962 -0.102623821
## Poatriv -0.1321800496 -0.45095346 0.103574859 -0.271062401 0.118784798
## Ranuflam 0.2792578969 0.03027933 -0.141161681 -0.103865168 0.001850673
## Rumeacet -0.1082864119 -0.02990082 0.094139189 -0.241619530 0.430903141
## Sagiproc 0.0410750255 -0.10193982 0.429563593 0.251953139 -0.030424463
## Salirepe 0.0654008835 0.19415423 0.012563845 0.269469132 -0.034009563
## Scorautu -0.0564104250 0.19693367 0.165511278 0.187848971 0.054307016
## Trifprat -0.1004765803 0.03522370 -0.042431386 -0.187403726 0.283990277
## Trifrepe -0.0472757003 -0.00868603 0.092022393 0.090353422 0.285088648
## Vicilath -0.0543676300 0.06549366 -0.040256932 0.226673409 0.078007690
## Bracruta 0.0855581430 0.13397561 -0.006591469 0.358036940 0.449131491
## Callcusp 0.2051900747 0.07088792 -0.181340715 -0.128624446 -0.029851304
## PC6 PC7 PC8 PC9 PC10
## Achimill -0.001433447 -4.245252e-05 -0.139275612 -0.03103964 -0.144243870
## Agrostol 0.046045557 -3.525837e-02 -0.078558371 -0.10461922 0.210297407
## Airaprae -0.089716769 -1.169179e-02 -0.018998063 -0.09251269 0.013279470
## Alopgeni -0.051261763 1.777523e-03 -0.380074914 -0.02631089 0.158579434
## Anthodor -0.112317801 -1.080043e-01 -0.106432326 -0.23023176 0.031200744
## Bellpere 0.213317439 -4.816241e-01 -0.271368034 -0.03193991 0.036080560
## Bromhord 0.223494039 -3.649116e-01 -0.185180260 -0.27427582 -0.154746474
## Chenalbu 0.041261458 6.203707e-02 -0.064291687 0.13439007 -0.078827421
## Cirsarve 0.015441982 -9.199161e-02 -0.008592696 -0.13102732 0.142525942
## Comapalu 0.253410826 9.179244e-02 0.230556545 -0.18290679 -0.056799466
## Eleopalu 0.008558863 1.001312e-01 -0.072559072 -0.25624816 -0.007112988
## Elymrepe -0.134141352 -3.202739e-01 0.473903145 -0.06448830 -0.019325678
## Empenigr -0.015937423 -1.006169e-01 0.116818638 -0.06552646 -0.010327426
## Hyporadi -0.077618969 1.398465e-01 0.031234356 -0.16995428 0.102773191
## Juncarti -0.373541740 1.366486e-02 -0.035830116 -0.26948199 -0.557137446
## Juncbufo -0.050382648 1.472281e-01 0.136186952 0.29199405 -0.302167976
## Lolipere -0.077527961 2.633129e-01 0.057940579 -0.13474471 0.104356918
## Planlanc -0.061099106 1.501356e-01 -0.257949680 0.08738485 0.137031409
## Poaprat -0.066342051 2.883555e-01 -0.173483287 0.08478524 -0.165446302
## Poatriv -0.032482359 -9.052647e-02 -0.169807161 0.02929114 -0.104717464
## Ranuflam 0.065379870 1.165436e-01 -0.097973268 0.12783117 -0.220685143
## Rumeacet -0.190843551 -8.174166e-02 0.357868885 -0.01191780 0.056890788
## Sagiproc -0.020650816 1.011335e-01 0.130867029 -0.25731110 0.034324439
## Salirepe -0.138332766 -3.903741e-01 0.014296262 0.48475738 -0.197219009
## Scorautu 0.245056856 -6.067166e-02 -0.046290901 -0.05730164 -0.333608160
## Trifprat -0.099328659 -1.856270e-02 0.144877989 -0.01182233 0.118607508
## Trifrepe 0.603878190 5.671622e-02 0.209187168 -0.13289241 -0.132042485
## Vicilath 0.081648423 1.144655e-01 -0.154486392 0.04787343 0.032832865
## Bracruta -0.332896166 -2.195302e-01 -0.081168769 -0.21669754 0.145243528
## Callcusp 0.137925615 -4.803121e-02 0.086693108 0.31331767 0.340411008
## PC11 PC12 PC13 PC14 PC15
## Achimill -1.255743e-02 0.046215103 -0.256723257 0.1313444344 -0.057272752
## Agrostol 2.532151e-02 -0.420175415 0.252885419 -0.0253551477 0.347506784
## Airaprae 9.391228e-02 0.060848890 0.077436256 -0.0404887435 0.106565696
## Alopgeni 2.913598e-01 0.340305846 -0.063431110 -0.2035060907 -0.108231449
## Anthodor -5.605106e-02 0.047468359 0.162324336 -0.0767398112 0.536035731
## Bellpere 2.277214e-01 -0.066340892 0.071266730 -0.0407483438 -0.002383024
## Bromhord -3.973954e-01 -0.221153049 -0.277170883 -0.0725695655 -0.125790738
## Chenalbu -5.683765e-02 -0.062470024 0.161943974 0.1147861266 -0.051639846
## Cirsarve -1.545945e-01 -0.284355335 0.102264776 -0.0094803775 0.017698853
## Comapalu 2.082232e-01 -0.122194920 0.081787344 0.3183693811 0.048304841
## Eleopalu -3.295834e-02 0.171174153 -0.031189575 -0.0927636652 -0.164919807
## Elymrepe 2.981940e-01 -0.141465524 -0.013771024 -0.1796833253 -0.059413586
## Empenigr -1.562428e-01 0.232654972 0.076458863 0.0736090906 0.104395306
## Hyporadi 6.753349e-05 0.072357556 -0.044834636 -0.1755997380 -0.171165046
## Juncarti 6.361967e-02 0.016314991 -0.103042130 -0.3378149280 0.160162220
## Juncbufo 2.254486e-02 -0.235107742 -0.357842523 0.0954648841 0.192810046
## Lolipere -2.665899e-01 0.181612509 0.033499240 -0.0317211176 0.151252550
## Planlanc 1.713019e-01 -0.346381680 0.025105687 -0.2449153380 -0.091817541
## Poaprat 8.882030e-02 -0.157511234 0.260442751 -0.2080073931 0.200952629
## Poatriv -2.214968e-01 0.194808307 0.181781411 0.0678421701 0.010998379
## Ranuflam -1.877691e-01 -0.031168655 0.298475438 0.0007910441 -0.153366894
## Rumeacet 7.013359e-02 -0.009512001 0.043474284 -0.2219769387 -0.145511466
## Sagiproc -4.214871e-01 -0.209789606 -0.013130855 -0.1631780935 -0.143417999
## Salirepe -2.200017e-01 0.113806094 0.160506778 -0.0424753661 0.127010328
## Scorautu 1.534068e-01 -0.042566093 0.244233847 -0.2416944201 -0.275228477
## Trifprat -1.120780e-01 0.073288375 0.385957343 -0.0618061194 -0.193278565
## Trifrepe 7.257752e-02 0.289166907 0.008895578 -0.1890818479 0.325637771
## Vicilath 1.591913e-02 -0.091882888 -0.262687703 -0.0417491662 0.062814819
## Bracruta 6.889939e-02 0.075253778 -0.096255783 0.2311884646 0.117003891
## Callcusp -1.820802e-01 0.060639205 -0.233039917 -0.5210466770 0.143774250
## PC16 PC17 PC18 PC19 PC20
## Achimill -0.267107657 -0.41210202 0.22489767 0.09626702 -0.07897943
## Agrostol -0.199229193 0.04314512 0.29537634 -0.10446748 0.06311722
## Airaprae -0.172178271 0.10130470 -0.19670085 0.19247304 0.18122527
## Alopgeni -0.192503646 0.01423069 -0.03731558 -0.11972890 0.14908814
## Anthodor 0.124188453 -0.05692455 -0.01459231 -0.33620497 -0.11244598
## Bellpere 0.118498260 -0.06025397 -0.15845467 -0.02686258 -0.42878365
## Bromhord -0.108289301 0.33193794 0.02391451 -0.01840574 0.14495526
## Chenalbu 0.287726031 0.05372402 0.11981393 0.20287560 -0.02388393
## Cirsarve -0.176766183 -0.07102134 -0.11240559 0.40560643 0.32331465
## Comapalu -0.047010757 0.12643349 -0.15710049 -0.11232080 -0.19504909
## Eleopalu 0.036796973 0.11125275 -0.49725476 -0.12181591 -0.03419597
## Elymrepe 0.255508256 0.03734245 0.05576828 -0.11750010 0.15496839
## Empenigr 0.143748846 0.27864541 -0.18233754 0.05932814 0.20581585
## Hyporadi 0.140080500 0.35880767 0.33217606 0.02822715 -0.13068216
## Juncarti 0.052524852 -0.06945109 0.11264250 0.22398356 -0.02128059
## Juncbufo -0.114034976 0.36869437 -0.13618729 -0.08144971 -0.16453618
## Lolipere -0.272395734 0.13416285 0.10854839 -0.33124013 0.02818100
## Planlanc 0.174097492 0.08161027 -0.18223409 -0.10838544 0.31479145
## Poaprat -0.036402138 0.05182841 -0.23735935 0.22250549 -0.20920638
## Poatriv 0.295850357 0.12398510 0.01209593 0.03105173 -0.03980571
## Ranuflam 0.217490753 -0.05986007 0.20438397 -0.27663575 0.14252283
## Rumeacet -0.104640190 -0.15313263 -0.04734271 -0.17877450 0.07037098
## Sagiproc 0.115654700 -0.39453594 -0.23802349 -0.10195268 -0.24920446
## Salirepe -0.161050639 -0.07780178 -0.13606538 -0.12265427 0.13147778
## Scorautu -0.280606173 0.13110453 0.18988616 -0.18533745 -0.08861247
## Trifprat -0.086183943 0.15164056 0.10448329 0.32730883 -0.23536861
## Trifrepe 0.084685940 -0.12711462 0.04985696 0.19051210 0.20983191
## Vicilath 0.392993308 -0.09705062 0.14360459 -0.04214081 0.11269105
## Bracruta 0.058168027 0.09580990 0.15422886 0.09217431 -0.10719603
## Callcusp -0.006585697 0.09775620 0.06929603 0.11123932 -0.28379314Ouverture des librairies nécessaires à la création du graphique. Si problème avec installation de la librairie ggord : https://rdrr.io/github/fawda123/ggord/
library(ggord)
library(ggplot2)
library(ggrepel) # pour éviter overlap du texteOn va d’abord commencer le graphique avec la fonctions ggord du package du même nom.
Le graphique créé par la fonction ggord est un objet ggplot. Cet objet (ici appelé pca.plot) peut être personnalisé à même la fonction ggord, et peut aussi servir de base sur laquelle ajouter des couches avec ggplot.
pca.plot <- ggord(pca.dune,
xlims=c(-1,1), ylims=c(-1,1)) # xlims et ylims définissent les limites dans l'affichage des axe (les valeurs entre parenthèses correspondent à l'échelle numérique des axes). On ajuste au fur et à mesure, au besoin.
pca.plot # À la fin de chaque chunk, l'objet graphique est appelé pour qu'on visualise le rendu à chaque fois.
On remarque que la fonctionggord inclut par défaut la proportion de la variance expliquée par les axes. C’est super!
Par défaut, la fonction ggord va mettre en graphique les deux premiers axes de l’ordination. Mais si on s’intéresse à différents axes, on peut spécifier lesquels mettre en graphique.
pca.plot <- ggord(pca.dune,
xlims=c(-1,1), ylims=c(-1,1),
axes=c(1,3)) # Ici, on met en graphique les axes 1 et 3.
pca.plot
On revient aux deux premiers axes. Faisons quelques petits changements. Les changements sont marqués d’un “#” à côté de la ligne.
pca.plot <- ggord(pca.dune,
xlims=c(-1.2,1.4), ylims=c(-0.75,0.9), # ajustement des limites
axes=c(1,2), # retour aux axes 1 et 2
size=3, # diminution de la taille des points (sites)
arrow=NULL,# enlever les flèches (qui viennent par défaut) pour les espèces
labcol="forestgreen") # changer la couleur du texte pour les espèces
pca.plot
Faisons d’autres changements!
pca.plot <- ggord(pca.dune,
xlims=c(-1.2,1.4), ylims=c(-0.75,0.9),
axes=c(1,2),
size=3,
obslab=TRUE, # pour passer de points pour les sites (par défaut) à texte (nom des sites)
arrow=NULL,
labcol="forestgreen")
pca.plot
Il y a beaucoup d’espèces affichées au milieu et c’est illisible. En plus, les espèces au milieu ne contribuent pas beaucoup à la distinction compositionnelle des sites. On va donc garder seulement les quelques espèces dont les scores ont la valeur absolue la plus élevée, pour chacun des deux axes, qu’on met en graphique (ici PC1 et PC2).
D’abord, on va créer un nouveau dataframe contenant seulement ces espèces, en utilisant quelques fonctions de la librairie dplyr.
library(dplyr)
A <- top_n(sp.scores.df, 3, PC1) # sélectionne les 3 espèces dont le score est le plus élevé sur l'axe 1.
B <- top_n(sp.scores.df, -3, PC1) # sélectionne les 3 espèces dont le score est le plus bas sur l'axe 1.
C <- top_n(sp.scores.df, 3, PC2) # sélectionne les 3 espèces dont le score est le plus élevé sur l'axe 2.
D <- top_n(sp.scores.df, -3, PC2) # sélectionne les 3 espèces dont le score est le plus bas sur l'axe 2.
sp.scores.skim <- bind_rows(A,B,C,D) # merge les sub dataframe A, B, C et D.
(sp.scores.skim <- distinct(sp.scores.skim)) # on ne garde que les rangées (espèces) qui sont uniques (pour éviter qu'une même espèce s'y retrouve en copies)## PC1 PC2 PC3 PC4 PC5
## Agrostol 0.41561895 -0.22347029 0.03498520 -0.040936713 0.005936338
## Eleopalu 0.38086726 0.08575416 -0.32385581 -0.126443878 0.029551309
## Ranuflam 0.27925790 0.03027933 -0.14116168 -0.103865168 0.001850673
## Lolipere -0.33774272 -0.21557574 -0.36010199 0.304110194 0.018829208
## Planlanc -0.25558044 0.24517118 -0.06008393 -0.078852590 0.300164812
## Poaprat -0.28182916 -0.17723703 -0.14368835 0.167481962 -0.102623821
## Anthodor -0.19694223 0.29608445 0.16687027 -0.353895610 -0.077737618
## Hyporadi -0.05609736 0.24551212 0.24231892 0.072606975 -0.291583945
## Alopgeni 0.13911504 -0.39519834 0.32482657 0.002621325 0.006164517
## Elymrepe -0.14616936 -0.30816535 -0.16422486 0.020107591 -0.302907981
## Poatriv -0.13218005 -0.45095346 0.10357486 -0.271062401 0.118784798
## PC6 PC7 PC8 PC9 PC10
## Agrostol 0.046045557 -0.035258372 -0.07855837 -0.10461922 0.210297407
## Eleopalu 0.008558863 0.100131193 -0.07255907 -0.25624816 -0.007112988
## Ranuflam 0.065379870 0.116543559 -0.09797327 0.12783117 -0.220685143
## Lolipere -0.077527961 0.263312916 0.05794058 -0.13474471 0.104356918
## Planlanc -0.061099106 0.150135587 -0.25794968 0.08738485 0.137031409
## Poaprat -0.066342051 0.288355460 -0.17348329 0.08478524 -0.165446302
## Anthodor -0.112317801 -0.108004332 -0.10643233 -0.23023176 0.031200744
## Hyporadi -0.077618969 0.139846549 0.03123436 -0.16995428 0.102773191
## Alopgeni -0.051261763 0.001777523 -0.38007491 -0.02631089 0.158579434
## Elymrepe -0.134141352 -0.320273944 0.47390315 -0.06448830 -0.019325678
## Poatriv -0.032482359 -0.090526465 -0.16980716 0.02929114 -0.104717464
## PC11 PC12 PC13 PC14 PC15
## Agrostol 2.532151e-02 -0.42017542 0.25288542 -0.0253551477 0.34750678
## Eleopalu -3.295834e-02 0.17117415 -0.03118957 -0.0927636652 -0.16491981
## Ranuflam -1.877691e-01 -0.03116866 0.29847544 0.0007910441 -0.15336689
## Lolipere -2.665899e-01 0.18161251 0.03349924 -0.0317211176 0.15125255
## Planlanc 1.713019e-01 -0.34638168 0.02510569 -0.2449153380 -0.09181754
## Poaprat 8.882030e-02 -0.15751123 0.26044275 -0.2080073931 0.20095263
## Anthodor -5.605106e-02 0.04746836 0.16232434 -0.0767398112 0.53603573
## Hyporadi 6.753349e-05 0.07235756 -0.04483464 -0.1755997380 -0.17116505
## Alopgeni 2.913598e-01 0.34030585 -0.06343111 -0.2035060907 -0.10823145
## Elymrepe 2.981940e-01 -0.14146552 -0.01377102 -0.1796833253 -0.05941359
## Poatriv -2.214968e-01 0.19480831 0.18178141 0.0678421701 0.01099838
## PC16 PC17 PC18 PC19 PC20
## Agrostol -0.19922919 0.04314512 0.29537634 -0.10446748 0.06311722
## Eleopalu 0.03679697 0.11125275 -0.49725476 -0.12181591 -0.03419597
## Ranuflam 0.21749075 -0.05986007 0.20438397 -0.27663575 0.14252283
## Lolipere -0.27239573 0.13416285 0.10854839 -0.33124013 0.02818100
## Planlanc 0.17409749 0.08161027 -0.18223409 -0.10838544 0.31479145
## Poaprat -0.03640214 0.05182841 -0.23735935 0.22250549 -0.20920638
## Anthodor 0.12418845 -0.05692455 -0.01459231 -0.33620497 -0.11244598
## Hyporadi 0.14008050 0.35880767 0.33217606 0.02822715 -0.13068216
## Alopgeni -0.19250365 0.01423069 -0.03731558 -0.11972890 0.14908814
## Elymrepe 0.25550826 0.03734245 0.05576828 -0.11750010 0.15496839
## Poatriv 0.29585036 0.12398510 0.01209593 0.03105173 -0.03980571Mise en graphique de ces espèces seulement
pca.plot <- ggord(pca.dune,
xlims=c(-1.2,1.4), ylims=c(-0.75,0.9),
axes=c(1,2),
size=3,
obslab=TRUE,
arrow=NULL,
txt=NULL, # enlever le texte (par défaut) pour les espèces
addpts=sp.scores.skim, # Ajouter les espèces aux plus hauts et faibles scores
addsize=3, # taille des points des espèces
addcol="forestgreen") # couleur des espèces
# on enlève l'argument labcol
pca.plot
pca.plot <- ggord(pca.dune,
xlims=c(-1.2,1.4), ylims=c(-0.75,0.9),
axes=c(1,2),
size=3,
obslab=TRUE,
arrow=NULL,
txt=NULL,
addpts=sp.scores.skim,
addsize=4, # augmenter la taille du texte des espèces
ptslab=TRUE, # changer de points à texte pour les espèces
addcol="forestgreen")
pca.plot
On va ajouter les résultats d’un envfit. D’abord, on fait le test, en s’assurant qu’on calcule la régression sur les deux axes qu’on a choisit de mettre en graphique (par défaut, 1 et 2).
(pca.envfit <- envfit(pca.dune, dune.env))##
## ***VECTORS
##
## PC1 PC2 r2 Pr(>r)
## A1 0.97222 0.23409 0.3622 0.015 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
##
## ***FACTORS:
##
## Centroids:
## PC1 PC2
## Moisture1 -0.3635 0.0257
## Moisture2 -0.2224 -0.0314
## Moisture4 0.1408 -0.3321
## Moisture5 0.4504 0.0872
## ManagementBF -0.3597 0.0173
## ManagementHF -0.1930 -0.0745
## ManagementNM 0.2330 0.3751
## ManagementSF 0.1077 -0.3217
## UseHayfield -0.0994 0.2242
## UseHaypastu -0.0203 -0.2180
## UsePasture 0.1716 0.0350
## Manure0 0.2330 0.3751
## Manure1 -0.2294 -0.0161
## Manure2 -0.2385 -0.0895
## Manure3 0.1969 -0.1644
## Manure4 -0.1812 -0.3955
##
## Goodness of fit:
## r2 Pr(>r)
## Moisture 0.5366 0.001 ***
## Management 0.4614 0.002 **
## Use 0.1794 0.134
## Manure 0.4531 0.012 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999Pour incorporer les vecteurs résultant de ce test, il faut d’abord les extraire et les mettre dans un dataframe.
vect <- data.frame(pca.envfit[["vectors"]][["arrows"]])
vect$variable.env<-rownames(vect) # il mettre le rowname(la variable environnementale en colonne)
vect## PC1 PC2 variable.env
## A1 0.9722159 0.2340862 A1pca.plot <- ggord(pca.dune,
xlims=c(-1.2,1.4), ylims=c(-0.75,0.9),
axes=c(1,2),
size=3,
obslab=TRUE,
arrow=NULL,
txt=NULL,
addpts=sp.scores.skim,
addsize=4,
ptslab=TRUE,
addcol="forestgreen")
pca.plot +
coord_fixed() + # ajout des vecteurs de envfit
geom_segment(data = vect,
aes(x = 0, xend = PC1, y = 0, yend = PC2),
arrow = arrow(length = unit(0.5, "cm")),
colour = "darkgrey") +
geom_text_repel(data = vect,
aes(x = PC1, y = PC2, label = variable.env, size = 3),
show.legend = FALSE)## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.
Maintenant, on va ajouter des ellipses autour des sites appartenant à des groupes similaires. La fonction ggord est à même de calculer ces ellipses avec intervalle de confiance spécifiée. On va aussi revenir à des symboles au lieu du nom des sites et changer la couleur du texte des espèces.
pca.plot <- ggord(pca.dune,
xlims=c(-1.2,1.4), ylims=c(-0.75,0.9),
axes=c(1,2),
grp_in=dune.env$Management, # ajout des ellipses en fonction du type de Management
grp_title = "Management",
ellipse_pro=0.95, # valeur de confiance pour les ellipses
alpha_el=0.3, # transparence des ellipses
size=3,
obslab=FALSE, # changement de texte à points pour sites
arrow=NULL,
txt=NULL,
addpts=sp.scores.skim,
addsize=4,
ptslab=TRUE,
addcol="black", # couleur des espèces
alpha = 1)
pca.plot +
scale_shape_manual('Groups', values = c(15,16,17,18)) + # Symboles pour les sites en fonction de leur groupe. Assurer vous d'avoir un nombre de symboles correspondant au nombre de groupes.
coord_fixed() +
geom_segment(data = vect,
aes(x = 0, xend = PC1, y = 0, yend = PC2),
arrow = arrow(length = unit(0.5, "cm")),
colour = "darkgrey") +
geom_text_repel(data = vect,
aes(x = PC1, y = PC2, label = variable.env, size = 3),
show.legend = FALSE)## Scale for shape is already present.
## Adding another scale for shape, which will replace the existing scale.
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.
Voici une légende indiquant quels symboles sont associées aux différentes valeurs numériques de la couche scale_shape_manual.
Changer les couleurs des groupes (types d’aménagement). On va extraire les codes des couleurs d’une palette de RColorBrewer avec la fonction brewer.pal
library(RColorBrewer)
pca.colors <- brewer.pal(n = 4, name = "Dark2") # Codes de couleur héxadécimauxOn va incorporer ces codes de couleurs dans un argument. On déplace la légende vers le haut.
pca.plot <- ggord(pca.dune,
xlims=c(-1.2,1.4), ylims=c(-0.75,0.9),
axes=c(1,2),
grp_in=dune.env$Management,
grp_title = "Management",
cols=pca.colors, # ajout des couleurs choisies
ellipse_pro=0.95,
alpha_el=0.3,
size=3,
obslab=FALSE,
arrow=NULL,
txt=NULL,
addpts=sp.scores.skim,
addsize=4,
ptslab=TRUE,
addcol="black",
alpha = 1)
pca.plot <- pca.plot + # ici j'enregistre l'ensemble des arguments précédents et suivants dans mon objet pca.plot
scale_shape_manual('Groups', values = c(15,16,17,18)) +
coord_fixed() +
geom_segment(data = vect,
aes(x = 0, xend = PC1, y = 0, yend = PC2),
arrow = arrow(length = unit(0.5, "cm")),
colour = "darkgrey") +
geom_text_repel(data = vect,
aes(x = PC1, y = PC2, label = variable.env, size = 3),
show.legend = FALSE) +
theme(legend.position = 'top') # mettre la légende en haut## Scale for shape is already present.
## Adding another scale for shape, which will replace the existing scale.
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.pca.plot
Créer un thème custom pour ggplot. Allez explorer les options (liste des arguments dans le fichier help, lorsque vous rentrez ?theme dans la console de R).
super_theme <- theme(
panel.background = element_blank(),
panel.grid = element_blank(), # grillage, ici aucun, mais on aurait pu mettre par ex. element_line("black"),
axis.line = element_line("black"),
text = element_text(size = 12),
axis.text = element_text(size = 10, colour = "gray25"), # taille des valeurs des axes
axis.title = element_text(size = 14, colour = "gray25"), # taille des titres des axes
legend.title = element_text(size = 14),
legend.text = element_text(size = 14),
legend.key = element_blank())Ajouter ce thème au graphique.
pca.plot + super_theme
On va rajouter du texte!
pca.plot +
annotate("text", x = 1.1, y = 0.7, label = "PERMANOVA\nR2=30%\n p=0.007") +
super_theme
Sauvegarder
Satisfait? On sauvegarde l’image en png. Par défaut, ggplot enregistre le dernier graphique créé. Ici on copie-colle le code du graphique ci haut, puis sauvegarde ce graphique avec la fonction ggsave dans un dossier qu’on a créé (figures/).
pca.plot +
annotate("text", x = 1.1, y = 0.7, label = "PERMANOVA\nR2=30%\n p=0.007") +
super_theme
ggsave("pca.png", path="./figures/")## Saving 7 x 5 in imageCet atelier n’est pas fait pour vous montrer à utiliser une esthétique particulière, mais bien pour vous familiariser avec les options (presque infinies…) permettant de personnaliser un graphique d’ordination. Allez voir des figures d’ordination dans la littérature. Essayez de reproduire l’esthétique qui vous semble la plus appropriée… ou, tout simplement, amusez-vous les artistes!!
CA
Préparation des données
Faire une CA avec la fonction cca.
ca.dune <- cca(dune)
ca.summary <- summary(ca.dune)
summary(ca.dune)##
## Call:
## cca(X = dune)
##
## Partitioning of scaled Chi-square:
## Inertia Proportion
## Total 2.115 1
## Unconstrained 2.115 1
##
## Eigenvalues, and their contribution to the scaled Chi-square
##
## Importance of components:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7
## Eigenvalue 0.5360 0.4001 0.2598 0.17598 0.14476 0.10791 0.09247
## Proportion Explained 0.2534 0.1892 0.1228 0.08319 0.06844 0.05102 0.04372
## Cumulative Proportion 0.2534 0.4426 0.5654 0.64858 0.71702 0.76804 0.81175
## CA8 CA9 CA10 CA11 CA12 CA13 CA14
## Eigenvalue 0.08091 0.07332 0.05630 0.04826 0.04125 0.03523 0.020529
## Proportion Explained 0.03825 0.03466 0.02661 0.02282 0.01950 0.01665 0.009705
## Cumulative Proportion 0.85000 0.88467 0.91128 0.93410 0.95360 0.97025 0.979955
## CA15 CA16 CA17 CA18 CA19
## Eigenvalue 0.014911 0.009074 0.007938 0.007002 0.003477
## Proportion Explained 0.007049 0.004290 0.003753 0.003310 0.001644
## Cumulative Proportion 0.987004 0.991293 0.995046 0.998356 1.000000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
##
##
## Species scores
##
## CA1 CA2 CA3 CA4 CA5 CA6
## Achimill -0.90859 0.08461 -0.58636 -0.008919 -0.660183 -0.18877
## Agrostol 0.93378 -0.20651 0.28165 0.024293 -0.139326 -0.02256
## Airaprae -1.00434 3.06749 1.33773 0.194305 -1.081813 -0.53699
## Alopgeni 0.40088 -0.61839 0.85013 0.346740 0.016509 0.10169
## Anthodor -0.96676 1.08361 -0.17188 0.459788 -0.607533 -0.30425
## Bellpere -0.50018 -0.35503 -0.15239 -0.704153 -0.058546 0.07308
## Bromhord -0.65762 -0.40634 -0.30685 -0.496751 -0.561358 0.07004
## Chenalbu 0.42445 -0.84402 1.59029 1.248755 -0.207480 0.87566
## Cirsarve -0.05647 -0.76398 0.91793 -1.175919 -0.384024 -0.13985
## Comapalu 1.91690 0.52150 -1.18215 -0.021738 -1.359988 1.31207
## Eleopalu 1.76383 0.34562 -0.57336 -0.002976 -0.332396 -0.14688
## Elymrepe -0.37074 -0.74148 0.26238 -0.566308 -0.270122 -0.72624
## Empenigr -0.69027 3.26420 1.95716 -0.176936 -0.073518 -0.16083
## Hyporadi -0.85408 2.52821 1.13951 -0.175115 -0.311874 0.11177
## Juncarti 1.27580 0.09963 -0.09320 0.005536 0.289410 -0.78247
## Juncbufo 0.08157 -0.68074 1.00545 1.078390 0.268360 0.24168
## Lolipere -0.50272 -0.35955 -0.21821 -0.474727 0.101494 -0.01594
## Planlanc -0.84058 0.24886 -0.78066 0.371149 0.271377 0.11989
## Poaprat -0.38919 -0.32999 -0.02015 -0.358371 0.079296 -0.05165
## Poatriv -0.18185 -0.53997 0.23388 0.178834 -0.155342 -0.07584
## Ranuflam 1.55886 0.30700 -0.29765 0.046974 -0.008747 -0.14744
## Rumeacet -0.65289 -0.25525 -0.59728 1.160164 0.255849 -0.32730
## Sagiproc 0.00364 0.01719 1.11570 0.066981 0.186654 0.32463
## Salirepe 0.61035 1.54868 0.04970 -0.607136 1.429729 -0.55183
## Scorautu -0.19566 0.38884 0.03975 -0.130392 0.141232 0.23717
## Trifprat -0.88116 -0.09792 -1.18172 1.282429 0.325706 -0.33388
## Trifrepe -0.07666 -0.02032 -0.20594 0.026462 -0.186748 0.53957
## Vicilath -0.61893 0.37140 -0.46057 -1.000375 1.162652 1.44971
## Bracruta 0.18222 0.26477 -0.16606 0.064009 0.576334 0.07741
## Callcusp 1.95199 0.56743 -0.85948 -0.098969 -0.556737 0.23282
##
##
## Site scores (weighted averages of species scores)
##
## CA1 CA2 CA3 CA4 CA5 CA6
## 1 -0.81167 -1.0826714 -0.14479 -2.10665 -0.39287 -1.83462
## 2 -0.63268 -0.6958357 -0.09708 -1.18695 -0.97686 0.06575
## 3 -0.10148 -0.9128732 0.68815 -0.68137 -0.08709 -0.28678
## 4 -0.05647 -0.7639784 0.91793 -1.17592 -0.38402 -0.13985
## 5 -0.95293 -0.1846015 -0.95609 0.86853 -0.34552 -0.98333
## 6 -0.85633 -0.0005408 -1.39735 1.59909 0.65494 -0.19386
## 7 -0.87149 -0.2547040 -0.86830 0.90468 0.17385 -0.03446
## 8 0.76268 -0.2968459 0.35648 -0.10772 0.17507 -0.36444
## 9 0.09693 -0.7864314 0.86492 0.40090 0.28704 -1.02783
## 10 -0.87885 -0.0353136 -0.82987 -0.68053 -0.75438 0.81070
## 11 -0.64223 0.4440332 -0.17371 -1.09684 1.37462 2.00626
## 12 0.28557 -0.6656161 1.64423 1.71496 0.65381 1.17376
## 13 0.42445 -0.8440195 1.59029 1.24876 -0.20748 0.87566
## 14 1.91996 0.5351062 -1.39863 -0.08575 -2.21317 2.43044
## 15 1.91384 0.5079036 -0.96567 0.04227 -0.50681 0.19370
## 16 2.00229 0.1090627 -0.33414 0.33760 -0.50097 -0.76159
## 17 -1.47545 2.7724102 0.40859 0.75117 -2.59425 -1.10122
## 18 -0.31241 0.6328355 -0.66501 -1.12728 2.65575 0.97565
## 19 -0.69027 3.2642026 1.95716 -0.17694 -0.07352 -0.16083
## 20 1.94438 1.0688809 -0.66595 -0.55317 1.59606 -1.70292Extraction de la proportion de la variance expliquée par les deux premiers axes (pour s’y référer plus rapidement).
(prop.expl.ca1 <- ca.summary$cont$importance[2, "CA1"])## [1] 0.2533987(prop.expl.ca2 <- ca.summary$cont$importance[2, "CA2"])## [1] 0.1891696envfit
(ca.envfit <- envfit(ca.dune, dune.env))##
## ***VECTORS
##
## CA1 CA2 r2 Pr(>r)
## A1 0.998160 0.060614 0.3104 0.043 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
##
## ***FACTORS:
##
## Centroids:
## CA1 CA2
## Moisture1 -0.7484 -0.1423
## Moisture2 -0.4652 -0.2156
## Moisture4 0.1827 -0.7315
## Moisture5 1.1143 0.5708
## ManagementBF -0.7258 -0.1413
## ManagementHF -0.3867 -0.2960
## ManagementNM 0.6517 1.4405
## ManagementSF 0.3376 -0.6761
## UseHayfield -0.2861 0.6488
## UseHaypastu -0.0735 -0.5602
## UsePasture 0.5163 0.0508
## Manure0 0.6517 1.4405
## Manure1 -0.4639 -0.1738
## Manure2 -0.5872 -0.3600
## Manure3 0.5187 -0.3172
## Manure4 -0.2059 -0.8775
##
## Goodness of fit:
## r2 Pr(>r)
## Moisture 0.4113 0.005 **
## Management 0.4441 0.002 **
## Use 0.1845 0.087 .
## Manure 0.4552 0.004 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999Pour incorporer les vecteurs résultant de ce test, il faut d’abord les extraire et les mettre dans un dataframe.
vect.ca.env <- data.frame(ca.envfit[["vectors"]][["arrows"]])
vect.ca.env$variable.env <- rownames(vect.ca.env)
vect.ca.env## CA1 CA2 variable.env
## A1 0.9981613 0.06061359 A1PERMANOVA. On va tester si la composition végétale diffère entre différents type d’aménagement (facteur; variable quantitative) des sites.
dune.dist.chi2 <- vegdist(dune, method='chisq')
disper.dune.chi2 <- betadisper(dune.dist.chi2, dune.env$Management)
anova(disper.dune.chi2)## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 3 3.0277 1.0092 5.9124 0.006493 **
## Residuals 16 2.7311 0.1707
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1permanova.dune.ca <- adonis2(dune.dist.chi2 ~ Management, data = dune.env)
permanova.dune.ca## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dune.dist.chi2 ~ Management, data = dune.env)
## Df SumOfSqs R2 F Pr(>F)
## Management 3 12.399 0.25508 1.8263 0.01 **
## Residual 16 36.208 0.74492
## Total 19 48.606 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ici, le résultat du test de dispersion indique que la différence dans la composition végétale entre les types d’aménagement qui est calculée par la PERMANOVA pourrait être du à la différence dans la variance interne en composition des groupes (les groupes ont une dispersion hétérogène).
Extraire les scores de CA avec la fonction scores. Pour spécifier à ggplot quels sont les “points” ou le “texte” à mettre en graphique.
site.scores.dune.ca <- scores(ca.dune, display="sites", choices=c(1,2)) # ici on garde que les scores des axes 1 et 2 mais on aurait pu extraire tous les scores aussi.
site.scores.dune.ca <- data.frame(site.scores.dune.ca)
Management <- dune.env$Management
(site.scores.ca <- cbind(site.scores.dune.ca, Management)) # incorporer variable Management dans les scores (pour aider plus loin dans ggplot à changer l'aspect visuel des points des sites en fonction de leur Management)## CA1 CA2 Management
## 1 -0.81167372 -1.0826713631 SF
## 2 -0.63267723 -0.6958357018 BF
## 3 -0.10147537 -0.9128731535 SF
## 4 -0.05647118 -0.7639783770 SF
## 5 -0.95292777 -0.1846015038 HF
## 6 -0.85632779 -0.0005407796 HF
## 7 -0.87149118 -0.2547039586 HF
## 8 0.76268386 -0.2968458686 HF
## 9 0.09693281 -0.7864314398 HF
## 10 -0.87885097 -0.0353135691 BF
## 11 -0.64223343 0.4440332251 BF
## 12 0.28557236 -0.6656161454 SF
## 13 0.42445102 -0.8440194551 SF
## 14 1.91996119 0.5351062139 NM
## 15 1.91384431 0.5079035784 NM
## 16 2.00229065 0.1090626959 SF
## 17 -1.47545057 2.7724102369 NM
## 18 -0.31241136 0.6328354520 NM
## 19 -0.69026877 3.2642026273 NM
## 20 1.94438046 1.0688808879 NMsp.scores.dune.ca <- scores(ca.dune, display="species", choices=c(1,2)) # idem scores axes 1 et 2 seulement
(sp.scores.dune.ca <- data.frame(sp.scores.dune.ca))## CA1 CA2
## Achimill -0.908593994 0.08460595
## Agrostol 0.933782632 -0.20651410
## Airaprae -1.004341487 3.06748567
## Alopgeni 0.400877498 -0.61838999
## Anthodor -0.966759907 1.08360766
## Bellpere -0.500177308 -0.35502842
## Bromhord -0.657624283 -0.40634288
## Chenalbu 0.424451025 -0.84401946
## Cirsarve -0.056471184 -0.76397838
## Comapalu 1.916902749 0.52150490
## Eleopalu 1.763825952 0.34562338
## Elymrepe -0.370742471 -0.74147804
## Empenigr -0.690268768 3.26420263
## Hyporadi -0.854079092 2.52821112
## Juncarti 1.275799635 0.09962851
## Juncbufo 0.081568571 -0.68074282
## Lolipere -0.502724848 -0.35955008
## Planlanc -0.840581620 0.24885595
## Poaprat -0.389190052 -0.32998545
## Poatriv -0.181852449 -0.53996706
## Ranuflam 1.558855994 0.30699556
## Rumeacet -0.652892711 -0.25524774
## Sagiproc 0.003639789 0.01718622
## Salirepe 0.610351084 1.54868352
## Scorautu -0.195656929 0.38883571
## Trifprat -0.881164096 -0.09792387
## Trifrepe -0.076660735 -0.02032460
## Vicilath -0.618932296 0.37139708
## Bracruta 0.182224012 0.26477423
## Callcusp 1.951985809 0.56742556Comme on a fait acvec les scores des espèces de la PCA, on ne garde que les plus hautes valeurs absolues.
A.ca <- top_n(sp.scores.dune.ca, 3, CA1)
B.ca <- top_n(sp.scores.dune.ca, -3, CA1)
C.ca <- top_n(sp.scores.dune.ca, 3, CA2)
D.ca <- top_n(sp.scores.dune.ca, -3, CA2)
sp.scores.skim.ca <- bind_rows(A.ca,B.ca,C.ca,D.ca)
(sp.scores.skim.ca <- distinct(sp.scores.skim.ca))## CA1 CA2
## Comapalu 1.91690275 0.52150490
## Eleopalu 1.76382595 0.34562338
## Callcusp 1.95198581 0.56742556
## Achimill -0.90859399 0.08460595
## Airaprae -1.00434149 3.06748567
## Anthodor -0.96675991 1.08360766
## Empenigr -0.69026877 3.26420263
## Hyporadi -0.85407909 2.52821112
## Chenalbu 0.42445102 -0.84401946
## Cirsarve -0.05647118 -0.76397838
## Elymrepe -0.37074247 -0.74147804Graphique simple avec les fonctions plot et ordiplot.
plot(ca.dune)
ordiellipse(ca.dune, dune.env$Management, display = "sites", kind = "sd", label = T)
ggplot
On commence le graphique ggplot. On peut faire un graphique qui ressemble à celui fait par plot en ajoutant “à la main” les sites et les espèces (en fonction de leurs scores extraits).
ca.plot <-
ggplot(data = site.scores.ca,
aes(x = CA1, y = CA2)) +
expand_limits(x=c(-4,4), y=c(-1,3)) + # sert à modifier l'échelle numérique du graphique
geom_text(data = site.scores.dune.ca,
aes(x = CA1, y = CA2, label= rownames(site.scores.dune.ca)),
color = "black") + # ajout du texte des sites
geom_text(data = sp.scores.dune.ca,
aes(x = CA1, y = CA2, label= rownames(sp.scores.dune.ca)),
color = "red") + # ajout du texte des espèces
xlab("CA1 (25.34%)") +
ylab("CA2 (18.92%)") # contrairement avec la fonction ggord, on écrit la variation expliqué par les axes soi même
ca.plot
On peut rajouter les vecteurs de note envfit.
ca.plot <-
ggplot(data = site.scores.ca,
aes(x = CA1, y = CA2)) +
expand_limits(x=c(-2,2.5), y=c(-1,3)) + # on zoom in sur l'échelle
geom_text(data = site.scores.dune.ca,
aes(x = CA1, y = CA2, label= rownames(site.scores.dune.ca)),
color = "black") +
geom_text(data = sp.scores.dune.ca,
aes(x = CA1, y = CA2, label= rownames(sp.scores.dune.ca)),
color = "red") +
xlab("CA1 (25.34%)") +
ylab("CA2 (18.92%)") +
geom_segment(data = vect.ca.env, # ajout des vecteurs envfit
aes(x = 0, xend = CA1, y = 0, yend = CA2),
arrow = arrow(length = unit(0.75, "cm")),
color = "black") +
geom_text_repel(data = vect.ca.env,
aes(x = CA1, y = CA2, label = variable.env, size = 3),
color = "firebrick1",
fontface = "bold",
show.legend = FALSE)
ca.plot
On va mettre des points pour les sites et changer leurs couleurs et leurs symboles en fonction de leur type d’aménagement. On va aussi mettre en graphique uniquement les espèces aux plus grands scores.
ca.plot <-
ggplot(data = site.scores.ca,
aes(x = CA1, y = CA2),
color = Management) + # on ajoute ici arg "color"
expand_limits(x=c(-2,2.5), y=c(-1,3)) +
geom_point(aes(shape = Management, color = Management),
size=3) + # on change de geom_text à geom_point pour les sites
xlab("CA1 (25.34%)") +
ylab("CA2 (18.92%)") +
geom_text_repel(data = sp.scores.skim.ca,
aes(x = CA1, y = CA2),
label=rownames(sp.scores.skim.ca),
color = "black") + # on garde geom_text pour les espèces mais on change de dataframe
geom_segment(data = vect.ca.env,
aes(x = 0, xend = CA1, y = 0, yend = CA2),
arrow = arrow(length = unit(0.75, "cm")),
color = "black") +
geom_text_repel(data = vect.ca.env,
aes(x = CA1, y = CA2, label = variable.env, size = 3),
color = "firebrick1",
fontface = "bold",
show.legend = FALSE) +
scale_shape_manual(values = c(15,16,17,18)) # on décide des symboles à utiliser (se référer à la charte plus haut)
ca.plot
On va ajouter des ellipses autour de nos sites en fonction de leur type d’aménagement. On va d’abord utiliser la fonction ggordiplot de la librairie ggordiplots pour extraire les données des ellipses. (on peut aussi extraire les données de hulls et spiders.) En même temps, ça nous donne un autre beau graphique, mais qui est plus complexe à modifier.
library(ggordiplots)## Loading required package: glueca.ggordiplot <- gg_ordiplot(ca.dune, groups = dune.env$Management, kind = "sd", conf = 0.95, pt.size = 3) # on spécifie l'intervalle de confiance des ellipses
names(ca.ggordiplot)## [1] "df_ord" "df_mean.ord" "df_ellipse" "df_hull" "df_spiders"
## [6] "plot"ca.ellipses <- ca.ggordiplot$df_ellipse # extraction des valeurs des ellipses
head(ca.ellipses)## Group x y
## 1 BF -0.5850575 0.7794705
## 2 BF -0.5691512 0.7777560
## 3 BF -0.5530958 0.7726193
## 4 BF -0.5369545 0.7640805
## 5 BF -0.5207911 0.7521735
## 6 BF -0.5046694 0.7369452On prend les valeurs des ellipses qu’on vient d’extraire pour mettre en graphique dans ggplot ces ellipses. Om change aussi les couleurs de nos groupes.
ca.plot <-
ggplot(data = site.scores.ca,
aes(x = CA1, y = CA2),
color = Management) +
expand_limits(x=c(-2,2.5), y=c(-1,3)) +
geom_point(aes(shape = Management, color = Management),
size=3) +
xlab("CA1 (25.34%)") +
ylab("CA2 (18.92%)") +
geom_text_repel(data = sp.scores.skim.ca,
aes(x = CA1, y = CA2),
label=rownames(sp.scores.skim.ca),
color = "black") +
geom_segment(data = vect.ca.env,
aes(x = 0, xend = CA1, y = 0, yend = CA2),
arrow = arrow(length = unit(0.75, "cm")),
color = "black") +
geom_text_repel(data = vect.ca.env,
aes(x = CA1, y = CA2, label = variable.env, size = 3),
color = "firebrick1",
fontface = "bold",
show.legend = FALSE) +
scale_shape_manual(values = c(15,16,17,18)) +
geom_path(data = ca.ellipses,
aes(x = x, y = y, color = Group),
show.legend = FALSE) + # ajout des ellipses
scale_color_brewer(palette="Dark2") # changement des couleurs des ellipses et des sites
ca.plot
Enfin, on pourrait rajouter notre thème personnalisé, mais optons pour un existant comme “bw”.
ca.plot + theme_bw()
Comme mentionnée précédemment, vous pouvez modifier presque à l’infini votre graphique!
On sauvegarde notre figure de CA.
ca.plot + theme_bw()
ggsave("ca.png", path="./figures/")## Saving 7 x 5 in imagePCoA
Préparation des données
Faire une PCoA avec la fonction prcomp.
dune.bray <- vegdist(dune, method='bray')Faire une PCoA dans R
pcoa.dune <- cmdscale(dune.bray, k =(nrow(dune) - 1), eig = TRUE)## Warning in cmdscale(dune.bray, k = (nrow(dune) - 1), eig = TRUE): only 14 of the
## first 19 eigenvalues are > 0summary(pcoa.dune)## Length Class Mode
## points 280 -none- numeric
## eig 20 -none- numeric
## x 0 -none- NULL
## ac 1 -none- numeric
## GOF 2 -none- numericDans cet exemple, on va s’intéresser aux axes 1 et 3, pour faire changement!
Extraire les proportions expliquées par les axes d’intérêt (1 et 3).
pcoa.eig2 <- eigenvals(pcoa.dune) # extraction des eigenvalues
pcoa.eig2 <- abs(pcoa.eig2) # calcul de leurs valeurs absolues
sum.pcoa.eig2 <- sum(pcoa.eig2) # total des eigenvalues
pcoa.prop.exp <- pcoa.eig2/sum.pcoa.eig2 # proportion expliquée par chacun des axes
(pcoa.axis1 <- pcoa.prop.exp[1])## [1] 0.3510557(pcoa.axis3 <- pcoa.prop.exp[3])## [1] 0.09439071Extraire les scores de la PCoA.
site.scores.pcoa <- scores(pcoa.dune) # site scores
sp.scores.pcoa <- wascores(pcoa.dune$points, dune) # species scores
site.scores.pcoa <- data.frame(site.scores.pcoa)
sp.scores.pcoa <- data.frame(sp.scores.pcoa)
Management <- dune.env$Management
site.scores.pcoa <- cbind(site.scores.pcoa, Management)Comme on a fait acec les scores des espèces de la PCA et de la CA, on ne garde que les plus hautes valeurs absolues, mais cette fois ci, des axes 1 et 3!.
A.pcoa <- top_n(sp.scores.pcoa, 3, X1)
B.pcoa <- top_n(sp.scores.pcoa, -3, X1)
C.pcoa<- top_n(sp.scores.pcoa, 3, X3)
D.pcoa <- top_n(sp.scores.pcoa, -3, X3)
sp.scores.skim.pcoa <- bind_rows(A.pcoa,B.pcoa,C.pcoa,D.pcoa)
(sp.scores.skim.pcoa <- distinct(sp.scores.skim.pcoa))## X1 X2 X3 X4 X5
## Comapalu 0.458911230 0.1005585008 0.125102902 -0.10064941 -0.116930204
## Eleopalu 0.439782223 -0.0007156372 0.106965962 -0.02841692 -0.008577587
## Callcusp 0.478165876 0.0484426585 0.130601445 -0.08277595 -0.033978450
## Achimill -0.288291335 0.0379745204 0.024675952 -0.09206847 -0.054125785
## Bromhord -0.257640216 -0.0863954310 0.007271944 -0.04584125 -0.050452440
## Trifprat -0.277616388 0.0548361248 -0.015862669 -0.04290624 -0.172247260
## Chenalbu 0.179872241 -0.2151706072 -0.273949689 -0.11812429 0.054398976
## Empenigr -0.005394285 0.4691160578 -0.217330420 0.14961117 0.161953074
## Juncbufo 0.059380874 -0.1672334185 -0.216754109 0.01068520 -0.006485779
## X6 X7 X8 X9 X10
## Comapalu -0.08308461 -0.034691754 -0.006325968 0.046184733 -0.002246251
## Eleopalu 0.03187448 0.003764263 0.018051782 -0.005171053 -0.011032580
## Callcusp -0.02214687 0.023556133 0.013738043 -0.002484335 0.032798306
## Achimill -0.00227727 0.014193770 0.015356868 0.001709240 -0.008461562
## Bromhord -0.05616890 0.070270108 0.040802596 0.013692138 -0.034231159
## Trifprat 0.15050556 -0.017699829 0.014358724 -0.002775301 0.025056357
## Chenalbu -0.09378649 -0.035811438 -0.138395894 -0.091709220 -0.037094509
## Empenigr -0.09775372 0.044962813 0.071770276 0.049827956 0.025898170
## Juncbufo 0.03276354 -0.054325318 -0.034914441 0.026597573 0.011229702
## X11 X12 X13 X14
## Comapalu 0.072444555 -0.006104436 0.0127702458 0.007221267
## Eleopalu -0.008371554 -0.014122056 0.0008705509 -0.002602251
## Callcusp -0.035297813 0.016298900 -0.0011275608 0.002066468
## Achimill -0.012001757 -0.001817249 -0.0039917245 -0.004955411
## Bromhord -0.007920231 0.005418645 0.0040940773 -0.003890661
## Trifprat -0.028308593 -0.011715116 0.0154123351 0.010475144
## Chenalbu -0.040870761 -0.002477992 0.0321341180 0.003153000
## Empenigr -0.041757420 -0.042319419 -0.0135562316 0.017897602
## Juncbufo -0.004033853 0.018627584 0.0021911465 -0.006478585envfit SUR LES AXES 1 ET 3 (on n’oublie pas!)
(pcoa.envfit <- envfit(pcoa.dune, dune.env, choices=c(1,3)))##
## ***VECTORS
##
## Dim1 Dim3 r2 Pr(>r)
## A1 0.9999700 0.0078085 0.379 0.019 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
##
## ***FACTORS:
##
## Centroids:
## Dim1 Dim3
## Moisture1 -0.2650 0.0607
## Moisture2 -0.1620 0.0000
## Moisture4 0.1065 -0.2369
## Moisture5 0.3271 0.0069
## ManagementBF -0.2694 0.0706
## ManagementHF -0.1350 -0.0394
## ManagementNM 0.1822 0.0370
## ManagementSF 0.0650 -0.0395
## UseHayfield -0.0608 -0.0240
## UseHaypastu -0.0227 -0.0210
## UsePasture 0.1213 0.0673
## Manure0 0.1822 0.0370
## Manure1 -0.1654 0.0175
## Manure2 -0.1648 -0.0944
## Manure3 0.1397 -0.0451
## Manure4 -0.1656 0.0944
##
## Goodness of fit:
## r2 Pr(>r)
## Moisture 0.6918 0.001 ***
## Management 0.2634 0.144
## Use 0.0614 0.690
## Manure 0.2892 0.211
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999Pour incorporer les vecteurs résultant de ce test, il faut d’abord les extraire et les mettre dans un dataframe.
vect.pcoa.env <- data.frame(pcoa.envfit[["vectors"]][["arrows"]])
vect.pcoa.env$variable.env <- rownames(vect.pcoa.env)
vect## PC1 PC2 variable.env
## A1 0.9722159 0.2340862 A1PERMANOVA. On va tester si la composition végétale diffère entre différents type d’aménagement (facteur; variable quantitative) des sites.
dune.dist.bray <- vegdist(dune, method='bray')
disper.dune.bray <- betadisper(dune.dist.bray, dune.env$Management)
anova(disper.dune.bray)## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 3 0.13831 0.046104 1.9506 0.1622
## Residuals 16 0.37816 0.023635permanova.dune.pcoa <- adonis2(dune.dist.bray ~ Management, data = dune.env)
permanova.dune.pcoa## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dune.dist.bray ~ Management, data = dune.env)
## Df SumOfSqs R2 F Pr(>F)
## Management 3 1.4686 0.34161 2.7672 0.001 ***
## Residual 16 2.8304 0.65839
## Total 19 4.2990 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Extraction des ellipses de la PCoA (axes 1 et 3!)
library(ggordiplots)
pcoa.ggordiplot <- gg_ordiplot(pcoa.dune, choices = c(1,3), groups = dune.env$Management, kind = "sd", conf = 0.95, pt.size = 3)
names(pcoa.ggordiplot)## [1] "df_ord" "df_mean.ord" "df_ellipse" "df_hull" "df_spiders"
## [6] "plot"pcoa.ellipses <- pcoa.ggordiplot$df_ellipseggplot
pcoa.plot <-
ggplot(data = site.scores.pcoa,
aes(x = Dim1, y = Dim3),
color = Management) +
expand_limits(x=c(-0.35,0.6), y=c(-0.35,0.35)) +
geom_point(aes(shape = Management, color = Management),
size=3) +
xlab("PCoA1 (35.10%)") +
ylab("PCoA3 (9.44%)") +
geom_text_repel(data = sp.scores.skim.pcoa,
aes(x = X1, y = X3),
label=rownames(sp.scores.skim.pcoa),
color = "black") +
geom_segment(data = vect.pcoa.env,
aes(x = 0, xend = Dim1, y = 0, yend = Dim3),
arrow = arrow(length = unit(0.75, "cm")),
color = "black") +
geom_text_repel(data = vect.pcoa.env,
aes(x = Dim1, y = Dim3, label = variable.env, size = 3),
color = "firebrick1",
fontface = "bold",
show.legend = FALSE) +
scale_shape_manual(values = c(15,16,17,18)) +
geom_path(data = pcoa.ellipses,
aes(x = x, y = y, color = Group),
show.legend = FALSE) +
scale_color_brewer(palette="Dark2") +
super_theme
pcoa.plot
ggsave("pcoa.png", path="./figures/")## Saving 7 x 5 in image